home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATHEMAT / 1769.ZIP / PFSAF.FOR < prev    next >
Text File  |  1988-10-05  |  54KB  |  2,750 lines
  1. $include: 'pfcom.for'
  2.     open(unit=7,file='algin',status='old')
  3.     open(unit=8,file='algout',status='unknown')
  4.     call datim
  5.     tbeg=tnow
  6.     do 4 i=1,nvar
  7. 4    ittb(i)=0
  8.     ios(1)=1
  9.     iosi(1)=1
  10.     ipv(1)=1
  11.     pv(1,1)='?'
  12.     ittb(1)=1
  13.     itte(1)=1
  14.     itb(1)=1
  15.     ite(1)=1
  16.     ivn(1)=2
  17.     ivp(1)=1
  18.     nb=7
  19.     nd=0
  20.     nv=1
  21.     ns=1
  22.     np=2
  23.     nt=1
  24.     maxe=9999
  25.     maxnt=0
  26.     maxite=0
  27.     itetop=1
  28.     ndump=0
  29.     minns=9999
  30.     lintot=0
  31.     idflg=1
  32.     call defh
  33.     ma=ntot
  34.     mb=nvar
  35.     mc=nterm
  36.     md=ninch
  37.     me=ninl
  38.     mg=nder
  39.     mh=nlab
  40.     mi=ndind
  41.     mj=mitay
  42.     write(8,1)ma,mb,mc,md,me,mg,mh,mi,mj
  43. 1    format(' version 85/6/21 ',
  44.      +   ' NTOT=',i6,' NVAR=',i4,' NTERM=',i5,' NINCH=',i5,/
  45.      +  ,' NINL=',i6,' NDER=',i3,' NLAB=',i3,' NDIND=',i2
  46.      +  ,' MITAY=',i3)
  47.     call monitr
  48.     call datim
  49.     write(8,90)nd,nv,ns,maxnt,maxite,lintot
  50. 90    format(' nu. derivatives=',i3,'  nu. variables=',i3,
  51.      +    '  nu. sums=',i3,/,' max nu. terms=',i5,
  52.      +  '  max total nu. variables=',i6,'  input lines=',i5)
  53.     t=tnow-tbeg
  54.     write(8,91)t
  55. 91    format(' time used=',f9.2,' seconds')
  56.     stop
  57.     end
  58.  
  59.  
  60.     subroutine add(is,js,t,ians)
  61. $include: 'pfcom.for'
  62.     nt1=nt+1
  63.     ib=ittb(is)
  64.     ie=itte(is)
  65.     jb=ittb(js)
  66.     je=itte(js)
  67.     i=itb(ie)
  68.     if(ivn(i).ne.2 .or. ivp(i).le.maxe)goto 3
  69.     do 1 it=ib,ie
  70.     i=itb(it)
  71.     if(ivn(i).eq.2 .and. ivp(i).gt.maxe)goto 2
  72. 1    continue
  73. 2    ie=it-1
  74. 3    i=itb(je)
  75.     if(ivn(i).ne.2 .or.ivp(i).le.maxe)goto 7
  76.     do 4 it=jb,je
  77.     i=itb(it)
  78.     if(ivn(i).eq.2 .and. ivp(i).gt.maxe)goto 5
  79. 4    continue
  80. 5    je=it-1
  81. 7    if(ib.gt.ie)goto 70
  82.     if(jb.gt.je)goto 12
  83.     if(icomp(ib,jb))8,60,30
  84. 8    jl=jb
  85.     if(jb.eq.je)goto 10
  86.     if(icomp(ib,je))10,16,18
  87. 10    call cpynt(jb,je,t)
  88. 12    call cpynt(ib,ie,1.d0)
  89. 14    call defsum(ians,nt1,nt)
  90.     return
  91. 16    save=tm(je)
  92.     tm(je)=tm(ib)/t+tm(je)
  93.     call cpynt(jb,je,t)
  94.     tm(je)=save
  95.     ib=ib+1
  96.     jb=je+1
  97.     goto 7
  98. 18    jr=je
  99. 20    if(jr-jl.eq.1)goto 28
  100.     jc=(jl+jr)/2
  101.     if(icomp(ib,jc))22,26,24
  102. 22    jl=jc
  103.     goto 20
  104. 24    jr=jc
  105.     goto 20
  106. 26    save=tm(jc)
  107.     tm(jc)=tm(ib)/t+tm(jc)
  108.     call cpynt(jb,jc,t)
  109.     tm(jc)=save
  110.     ib=ib+1
  111.     jb=jc+1
  112.     goto 7
  113. 28    call cpynt(jb,jl,t)
  114.     jb=jr
  115. 30    il=ib
  116.     if(ib.eq.ie)goto 31
  117.     if(icomp(jb,ie))31,36,38
  118. 31    call cpynt(ib,ie,1.d0)
  119. 32    call cpynt(jb,je,t)
  120.     goto 14
  121. 36    save=tm(ie)
  122.     tm(ie)=tm(ie)+t*tm(jb)
  123.     call cpynt(ib,ie,1.d0)
  124.     tm(ie)=save
  125.     jb=jb+1
  126.     ib=ie+1
  127.     goto 70
  128. 38    ir=ie
  129. 40    if(ir-il.eq.1)goto 48
  130.     ic=(il+ir)/2
  131.     if(icomp(jb,ic))42,46,44
  132. 42    il=ic
  133.     goto 40
  134. 44    ir=ic
  135.     goto 40
  136. 46    save=tm(ic)
  137.     tm(ic)=tm(ic)+t*tm(jb)
  138.     call cpynt(ib,ic,1.d0)
  139.     tm(ic)=save
  140.     jb=jb+1
  141.     ib=ic+1
  142.     goto 7
  143. 48    call cpynt(ib,il,1.d0)
  144.     ib=ir
  145.     goto 8
  146. 60    save=tm(ib)
  147.     tm(ib)=tm(ib)+t*tm(jb)
  148.     call cpynt(ib,ib,1.d0)
  149.     tm(ib)=save
  150.     ib=ib+1
  151.     jb=jb+1
  152.     goto 7
  153. 70    if(jb.gt.je)goto 14
  154.     goto 32
  155.     end
  156.  
  157.  
  158.     subroutine bomb(i,h)
  159. $include: 'pfcom.for'
  160.     character*8 h
  161.     goto(100,100,30,40,50),i
  162. 30    write(8,31)h
  163. 31    format(a8)
  164.     goto 100
  165. 40    write(8,41)h,h
  166. 41    format(a8,' is too small for your problem.',
  167.      +       '  Recompile with bigger ',a8)
  168.     goto 100
  169. 50    continue
  170. 100    write(8,101)line,ap
  171. 101    format('input line number',i5,' was:',/,127a1)
  172.     call datim
  173.     write(8,90)nd,nv,ns,maxnt,maxite,line
  174. 90    format(/,'nu. derivatives=',i3,
  175.      +      '   nu. variables=',i3,'nu. sums=',i3,'  max nu. terms=',
  176.      +      i5,/,'  max total nu. variables=',i6,' input lines=',i5)
  177.     if(ndump.ne.0)call dump
  178.     t=tnow-tbeg
  179.     write(8,91)t
  180. 91    format('time used=',f9.2,' seconds')
  181.     stop
  182.     end
  183.  
  184.  
  185.     subroutine bkup(j,k)
  186. $include: 'pfcom.for'
  187.     if(j.eq.k)goto 2
  188.     do 1 i=j,k-1
  189.     itb(i)=itb(i+1)
  190.     ite(i)=ite(i+1)
  191. 1    tm(i)=tm(i+1)
  192. 2    k=k-1
  193.     return
  194.     end
  195.  
  196.  
  197.     subroutine call
  198. $include: 'pfcom.for'
  199.     call nxt(3,nchpl,' ',j1)
  200.     call numl(j1,nchpl,j2,j3,j4)
  201.     if(linlab(j4).eq.-1)call bomb(3,'call    ')
  202.     linret=line
  203.     line=linlab(j4)
  204.     return
  205.     end
  206.  
  207.  
  208.     subroutine coef(is,j3,c)
  209. c  c=coefficient, j3=start of variable (on return)
  210. $include: 'pfcom.for'
  211.     f=1.d0
  212.     sign=1.d0
  213.     call nxtnb(is,nchpl,j2)
  214.     call srch(j2,nchpl,'*',jj)
  215.     if(jj.eq.-1)jj=nchpl
  216.     call srch(j2,nchpl,'/',jk)
  217.     if(jk.eq.-1)jk=nchpl
  218.     call srchch(j2,nchpl,j3)
  219.     j3=min0(jj,jk,j3)
  220.     if(a(j2).eq.'-')goto 4
  221.     if(a(j2).eq.'+')goto 10
  222.     j4=j2
  223.     goto 12
  224. 4    sign=-1.d0
  225. 10    call nxtnb(j2+1,nchpl,j4)
  226. 12    c=sign
  227.     if(a(j4).eq.'.')goto 15
  228.     call intval(j4,i)
  229.     if(i.eq.-1 .and. j4.ne.j3)call bomb(3,'coef1   ')
  230.     if(i.eq.-1)return
  231. 15    c=0.d0
  232.     do 30 j1=j4,nchpl
  233.     if(a(j1).eq.'.' .and. f.lt..5d0)call bomb(3,'coef2   ')
  234.     if(a(j1).eq.'.')f=.1d0
  235.     if(a(j1).eq.'.')goto 30
  236.     call intval(j1,i)
  237.     if(i.eq.-1)goto 40
  238.     if(f.gt..5d0)c=10.d0*c+i
  239.     if(f.lt..5d0)c=c+i*f
  240.     if(f.lt..5d0)f=f*.1d0
  241. 30    continue
  242. 40    c=c*sign
  243.     return
  244.     end
  245.  
  246.  
  247.     subroutine comput
  248. $include: 'pfcom.for'
  249. 300    call nxt(1,nchpl,' ',j1)
  250.     call numv(j1,nchpl,j2,j3,j4)
  251.     call nxt(j3,nchpl,'=',j5)
  252.     call srch(j5,nchpl,'/',j6)
  253.     if(j6.eq.-1)goto 320
  254. c  /
  255.     call nxtnb(j5+1,j6-1,j7)
  256.     if(a(1).eq.'d')goto 301
  257.     call numv(j7,j6-1,j8,j9,j10)
  258.     call numv(j6+1,nchpl,j11,j12,j13)
  259.     if(ittb(j10).eq.0 .or. ittb(j13).eq.0)call bomb(3,'not sum ')
  260.     call recip(j13,1)
  261.     nt1=nt+1
  262.     call mult(ittb(j10),itte(j10),ittb(1),itte(1))
  263.     call defsum(j4,nt1,nt)
  264.     call outqu(j4)
  265.     minns=min0(minns,ios(1))
  266.     it=ittb(1)
  267.     itte(1)=it
  268.     ite(it)=itb(it)
  269.     k=itb(it)
  270.     ivn(k)=2
  271.     ivp(k)=0
  272.     call pack
  273.     return
  274. c  compute sj4=dsj10/dfj14
  275. 301    call numv(j7+1,j6-1,j8,j9,j10)
  276.     if(a(j7).ne.'d')call bomb(3,'comput1 ')
  277.     if(idflg.ne.0)call dertab
  278.     call nxt(j6,nchpl,'d',j11)
  279.     call numv(j11+1,nchpl,j12,j13,j14)
  280.     call srch(j13,nchpl,'=',j15)
  281.     iord=1
  282.     if(j15.eq.-1)goto 311
  283.     call nxtnb(j15+1,nchpl,j16)
  284.     call intval(j16,j17)
  285.     if(j17.ne.-1)goto 360
  286.     call numv(j16,nchpl,j19,j20,j18)
  287.     if(ittb(j18).eq.0)call bomb(3,'comput14')
  288.     iord=1.000000001d0*tm(ittb(j18))
  289.     if(iord.gt.0)goto 311
  290.     call bomb(3,'comput15')
  291. 360    call getexp(j15+1,iord,j16,0)
  292.     if(j16.ne.-1 .or. iord.le.0)call bomb(3,'comput2 ')
  293. 311    do 313 ii=1,iord
  294.     if(ii.eq.1)call deriv(j10,j14,1)
  295.     if(ii.ne.1)call deriv(1,j14,1)
  296.     call outqu(1)
  297.     if(tm(ittb(1)).eq.0.)goto 314
  298. 313    continue
  299. 314    call swap(j4)
  300.     return
  301. c  not /
  302. 320    call srch(j5,nchpl,'+',j6)
  303.     if(a(1).eq.'d')call bomb(3,'comput3 ')
  304.     if(j6.ne.-1)goto 321
  305.     call srch(j5,nchpl,'-',j6)
  306.     tmm=-1.d0
  307.     if(j6.ne.-1)goto 322
  308.     call srch(j5,nchpl,'*',j6)
  309.     if(j6.ne.-1)goto 339
  310.     call srch(j5,nchpl,'^',j6)
  311.     if(j6.ne.-1)goto 339
  312. c  compute sj4=sj12
  313.     call numv(j5,nchpl,j6,j11,j12)
  314.     if(ittb(j12).eq.0)call bomb(3,'comput4 ')
  315.     call copys(j12,j4)
  316.     call outqu(j4)
  317.     return
  318. c  compute sj4=sj9  + or -  sj12
  319. 321    tmm=1.d0
  320. 322    call numv(j5+1,j6-1,j7,j8,j9)
  321.     call numv(j6+1,nchpl,j10,j11,j12)
  322.     call add(j9,j12,tmm,j4)
  323.     call outqu(j4)
  324.     return
  325. c  compute sj4=sj9*sj12
  326. 339    call numv(j5+1,j6-1,j7,j8,j9)
  327.     if(a(j6).eq.'^')goto 350
  328.     if(a(j6+1).eq.'*')goto 349
  329.     call numv(j6+1,nchpl,j10,j11,j12)
  330.     nt1=nt+1
  331.     call mult(ittb(j9),itte(j9),ittb(j12),itte(j12))
  332.     call defsum(j4,nt1,nt)
  333.     call outqu(j4)
  334.     return
  335. 349    j6=j6+1
  336. 350    call numv(j6+1,nchpl,j10,j11,j12)
  337.     call expnd(j9,j12,j4)
  338.     return
  339.     end
  340.  
  341.  
  342.     subroutine copy(a,b,n)
  343.     character*1 a,b
  344.     dimension a(2),b(2)
  345.     do 1 i=1,n
  346. 1    b(i)=a(i)
  347.     return
  348.     end
  349.  
  350.     subroutine copys(j1,j2)
  351. $include: 'pfcom.for'
  352.     nt1=nt+1
  353.     call cpynt(ittb(j1),itte(j1),1.d0)
  354.     call defsum(j2,nt1,nt)
  355.     return
  356.     end
  357.  
  358.  
  359.     subroutine count
  360. $include: 'pfcom.for'
  361.     call nxt(3,nchpl,' ',j1)
  362.     call numv(j1,nchpl,j3,j4,j5)
  363.     do 12 i=nchpl,j4+1,-1
  364.     if(a(i).ne.' ')goto 14
  365. 12    continue
  366.     call bomb(3,'count1  ')
  367. 14    do 16 j=i,j4+1,-1
  368.     if(a(j).eq.' ')goto 18
  369. 16    continue
  370. 18    call numv(j,i,j6,j7,j8)
  371.     if(ittb(j8).eq.0)call bomb(3,'count2  ')
  372.     if(ittb(j5).ne.0)goto 20
  373.     k=itetop+1
  374.     call incnt
  375.     itb(nt)=k
  376.     call nite(k)
  377.     tm(nt)=itte(j8)+1-ittb(j8)
  378.     ivn(k)=2
  379.     ivp(k)=0
  380.     call defsum(j5,nt,nt)
  381.     return
  382. 20    ia=ittb(j5)
  383.     ib=itte(j5)
  384.     k=itb(ia)
  385.     l=ite(ia)
  386.     ivn(k)=2
  387.     ivp(k)=0
  388.     itte(j5)=ia
  389.     ite(ia)=k
  390.     tm(ia)=itte(j8)+1-ittb(j8)
  391.     if(ia.ne.ib .or. k.ne.l)minns=min0(minns,ios(j5))
  392.     call pack
  393.     return
  394.     end
  395.  
  396.  
  397.     subroutine cpynt(ia,ib,tmm)
  398. c  copy terms ia-ib to nt+1,..., coefficients are multiplied by tmm
  399. $include: 'pfcom.for'
  400.     if(ib.lt.ia)return
  401.     if(ia.lt.1 .or. ib.gt.nt)call bomb(3,'cpynt2  ')
  402.     do 15 i=ia,ib
  403.     ja=itb(i)
  404.     if(ivn(ja).eq.2 .and. ivp(ja).gt.maxe)goto 15
  405.     k=itetop
  406.     call incnt
  407.     tm(nt)=tm(i)*tmm
  408.     itb(nt)=k+1
  409.     do 10 j=ja,ite(i)
  410.     k=k+1
  411.     ivn(k)=ivn(j)
  412. 10    ivp(k)=ivp(j)
  413.     call nite(k)
  414. 15    continue
  415. 16    return
  416.     end
  417.  
  418.  
  419.     subroutine datim
  420. $include: 'pfcom.for'
  421.     call date(iy,mon,id)
  422.     call time(ih,min,is,if)
  423.     tnow=is+60.*(min+60.*ih)+if/100.
  424.     write(8,1)iy,mon,id,ih,min,is
  425. 1    format(' job run  ',i4,'/',i2,'/',i2,i9,':',i2,':',i2)
  426.     return
  427.     end
  428.  
  429.  
  430.     subroutine datime
  431. $include: 'pfcom.for'
  432.     call time(ih,min,is,if)
  433.     tnow=is+60.*(min+60.*ih)+if/100.
  434.     return
  435.     end
  436.  
  437.  
  438.     subroutine defc(a,n)
  439.     character*1 a(2)
  440.     do 1 i=2,n
  441.     j=a(i-1)
  442. 1    a(i)=j+1
  443.     return
  444.     end
  445.  
  446.  
  447.     subroutine defh
  448. $include: 'pfcom.for'
  449.     character*1 ga(26),gc(26),gn(10)
  450.     data ga/'a','b','c','d','e','f','g','h','i','j','k','l','m'
  451.      ,       ,'n','o','p','q','r','s','t','u','v','w','x','y','z'/
  452.     data gc/'A','B','C','D','E','F','G','H','I','J','K','L','M'
  453.      ,       ,'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
  454.     data gn/'0','1','2','3','4','5','6','7','8','9'/
  455.     do 10 i=1,10
  456. 10    hn(i)=gn(i)
  457.     do 20 i=1,26
  458.     ha(i)=ga(i)
  459. 20    hc(i)=gc(i)
  460.     return
  461.     end
  462.  
  463.  
  464.     subroutine defsum(is,ittbv,ittev)
  465. $include: 'pfcom.for'
  466.     if(ittev.ne.nt .or. ittev+1.ne.ittbv)goto 1
  467.     call niltrm
  468.     ittev=nt
  469. 1    if(ittbv.gt.ittev)call bomb(3,'defsum1 ')
  470.     if(ittb(is).eq.0)goto 20
  471.     j=ios(is)
  472.     minns=min0(minns,j)
  473.     if(ittbv.ge.ittb(is) .and. ittev.le.itte(is))goto 8
  474.     if(j.eq.ns .and. ittbv.ge.ittb(is))goto 8
  475.     if(ittbv.gt.itte(iosi(ns)))goto 4
  476.     call bomb(3,'defsum2 ')
  477. 4    ios(is)=ns
  478.     do 6 i=j,ns-1
  479.     k=iosi(i+1)
  480.     iosi(i)=k
  481. 6    ios(k)=i
  482.     iosi(ns)=is
  483. 8    ittb(is)=ittbv
  484.     itte(is)=ittev
  485.     call pack
  486.     return
  487. 20    if(ittbv.le.itte(iosi(ns)))call bomb(3,'defsum3 ')
  488.     ns=ns+1
  489.     ios(is)=ns
  490.     iosi(ns)=is
  491.     ittb(is)=ittbv
  492.     itte(is)=ittev
  493.     minns=min0(minns,ns)
  494.     call pack
  495.     return
  496.     end
  497.  
  498.  
  499.     subroutine delet
  500. $include: 'pfcom.for'
  501. c  delete sj4
  502.     call nxt(4,nchpl,' ',j1)
  503. 1    call numv(j1,nchpl,j2,j3,j4)
  504.     call delsum(j4)
  505.     if(j3.eq.nchpl)goto 10
  506.     call srchnb(j3+1,nchpl,j1)
  507.     if(j1.ne.-1)goto 1
  508. 10    return
  509.     end
  510.  
  511.  
  512.     subroutine delsum(is)
  513. $include: 'pfcom.for'
  514.     if(ittb(is).eq.0)return
  515.     j=ios(is)
  516.     minns=min0(minns,j)
  517.     if(j.eq.0 .or. ns.eq.0)call bomb(3,'delsum1 ')
  518.     ns=ns-1
  519.     if(j.gt.ns)goto 10
  520.     do 2 i=j,ns
  521.     k=iosi(i+1)
  522.     iosi(i)=k
  523. 2    ios(k)=i
  524. 10    iosi(ns+1)=0
  525.     ios(is)=0
  526.     ittb(is)=0
  527.     itte(is)=0
  528.     call pack
  529.     return
  530.     end
  531.  
  532.  
  533.     subroutine deriv(j10,ind,jans)
  534. c  d(sum j10)/d(var ind) is put in sum j4
  535. c  ider(1,nd) is dependent variable in derivative
  536. c  ider(2,nd) is independent variable in derivative
  537. c  ider(3,nd) is variable equal to derivative
  538. $include: 'pfcom.for'
  539.     character*1 bl
  540.     common/i/if
  541.     nt1=nt+1
  542.     if(ittb(j10).lt.1 .or. itte(j10).gt.nt)call bomb(3,'deriv2  ')
  543.     do 90 it=ittb(j10),itte(j10)
  544.     ia=itb(it)
  545.     ib=ite(it)
  546.     do 85 ic=ia,ib
  547.     ie=ivn(ic)
  548.     if(ie.ne.ind)goto 40
  549.     call incnt
  550.     itb(nt)=itetop+1
  551.     tm(nt)=tm(it)*ivp(ic)
  552.     if=itetop
  553.     do 36 i=ia,ib
  554.     if(i.eq.ic)goto 32
  555.     if=if+1
  556.     ivn(if)=ivn(i)
  557.     ivp(if)=ivp(i)
  558.     goto 36
  559. 32    if(ivp(i).eq.1 .and. ia.ne.ib)goto 36
  560.     if=if+1
  561.     ivn(if)=ivn(i)
  562.     ivp(if)=ivp(i)-1
  563. 36    continue
  564.     call nite(if)
  565.     if(if.lt.itb(nt))call bomb(3,'deriv3  ')
  566.     goto 85
  567. 40    if(nd.eq.0)goto 85
  568.     i1=1
  569. 41    j1=idt(i1,ie)
  570.     if(j1.eq.0)goto 85
  571.     k1=ider(2,j1)
  572.     level=1
  573.     if(k1.eq.ind)goto 50
  574.     i2=1
  575. 42    j2=idt(i2,k1)
  576.     if(j2.eq.0)goto 77
  577.     k2=ider(2,j2)
  578.     level=2
  579.     if(k2.eq.ind)goto 50
  580.     i3=1
  581. 43    j3=idt(i3,k2)
  582.     if(j3.eq.0)goto 76
  583.     k3=ider(2,j3)
  584.     level=3
  585.     if(k3.eq.ind)goto 50
  586.     i4=1
  587. 44    j4=idt(i4,k3)
  588.     if(j4.eq.0)goto 75
  589.     k4=ider(2,j4)
  590.     level=4
  591.     if(k4.eq.ind)goto 50
  592.     i5=1
  593. 45    j5=idt(i5,k4)
  594.     if(j5.eq.0)goto 74
  595.     k5=ider(2,j5)
  596.     level=5
  597.     if(k5.eq.ind)goto 50
  598.     l1=ipv(ie)
  599.     l2=ipv(k1)
  600.     l3=ipv(k2)
  601.     l4=ipv(k3)
  602.     l5=ipv(k4)
  603.     l6=ipv(k5)
  604.     bl=' '
  605.     write(8,49)(pv(i,ie),i=1,l1),bl,(pv(i,k1),i=1,l2),bl
  606.      +  ,(pv(i,k2),i=1,l3),bl,(pv(i,k3),i=1,l4),bl,(pv(i,k4),i=1,l5)
  607.      +  ,bl,(pv(i,k5),i=1,l6)
  608. 49    format('chain too long: ',90a1)
  609.     call bomb(3,'deriv4  ')
  610. 50    call incnt
  611.     if=itetop
  612.     itb(nt)=if+1
  613.     tm(nt)=tm(it)*ivp(ic)
  614.     do 55 i=ia,ib
  615.     if(i.eq.ic)goto 52
  616.     if=if+1
  617.     ivn(if)=ivn(i)
  618.     ivp(if)=ivp(i)
  619.     goto 55
  620. 52    if(ivp(i).eq.1)goto 55
  621.     if=if+1
  622.     ivn(if)=ivn(i)
  623.     ivp(if)=ivp(i)-1
  624. 55    continue
  625.     call deriv2(j1)
  626.     if(level.eq.1)goto 77
  627.     call deriv2(j2)
  628.     if(level.eq.2)goto 76
  629.     call deriv2(j3)
  630.     if(level.eq.3)goto 75
  631.     call deriv2(j4)
  632.     if(level.eq.4)goto 74
  633.     call deriv2(j5)
  634.     i5=i5+1
  635.     if(i5.le.ndind)goto 45
  636. 74    i4=i4+1
  637.     if(i4.le.ndind)goto 44
  638. 75    i3=i3+1
  639.     if(i3.le.ndind)goto 43
  640. 76    i2=i2+1
  641.     if(i2.le.ndind)goto 42
  642. 77    i1=i1+1
  643.     if(i1.le.ndind)goto 41
  644. 85    continue
  645. 90    continue
  646.     call defsum(jans,nt1,nt)
  647.     call order(jans)
  648.     return
  649.     end
  650.  
  651.  
  652.     subroutine deriv2(j)
  653. $include: 'pfcom.for'
  654.     common/i/if
  655.     if=if+1
  656.     call nite(if)
  657.     ivn(if)=ider(3,j)
  658.     ivp(if)=1
  659.     return
  660.     end
  661.  
  662.  
  663.     subroutine dertab
  664. $include: 'pfcom.for'
  665.     idflg=0
  666.     do 302 i=1,ndind
  667.     do 302 j=1,nvar
  668. 302    idt(i,j)=0
  669.     do 304 i=1,nd
  670.     idepv=ider(1,i)
  671.     do 303 j=1,ndind
  672.     if(idt(j,idepv).eq.0)goto 304
  673. 303    continue
  674.     call bomb(4,'ndind   ')
  675. 304    idt(j,idepv)=i
  676.     return
  677.     end
  678.  
  679.  
  680.     subroutine dfine
  681. $include: 'pfcom.for'
  682.     call nxt(4,nchpl,' ',j1)
  683.     call nxtch(j1,nchpl,j2)
  684.     call numv(j2,nchpl,j4,j5,j6)
  685.     ittb6=nt+1
  686.     call srch(j5,nchpl,'=',j7)
  687.     if(j7.eq.-1)goto 410
  688.     call term(j7+1)
  689.     goto 412
  690. 410    call read(ib)
  691.     if(np.ge.2)call prntap
  692.     if(ib.eq.1)goto 412
  693.     call term(1)
  694.     goto 410
  695. 412    call defsum(1,ittb6,nt)
  696.     call order(1)
  697.     call outqu(1)
  698.     call swap(j6)
  699.     return
  700.     end
  701.  
  702.     subroutine dump
  703. $include: 'pfcom.for'
  704.     character*1 bl
  705.     data bl/' '/
  706.     write(8,1)nd,nv,ns,np,nt,line,maxe,maxnt,maxite
  707. 1    format(/,'nd=',i3,' nv=',i3,' ns=',i3,' np=',i2,' nt=',i5,
  708.      +     ' line=',i3,' maxe=',i2,' maxnt=',i5,' maxite=',i6)
  709.     do 10 i=1,nv
  710.     j=ipv(i)
  711.     if(j.eq.nch)goto 3
  712.     jp=j+1
  713.     do 4 jj=jp,nch
  714. 4    pv(jj,i)=bl
  715. 3    write(8,5)i,(pv(k,i),k=1,nch),ittb(i),itte(i),ios(i),iosi(i)
  716. 5    format('i=',i5,'   pv=',15a1,'   ittb=',i5,'   itte=',i5
  717.      +      ,'  ios=',i3,'  iosi=',i3)
  718. 10    continue
  719.     write(8,11)
  720. 11    format(' ')
  721.     write(8,13)
  722. 13    format('       i     tm     itb   ite')
  723.     do 20 i=1,nt,4
  724.     j=i+1
  725.     k=i+2
  726.     l=i+3
  727.     write(8,21)i,tm(i),itb(i),ite(i),j,tm(j),itb(j),ite(j)
  728.      +          ,k,tm(k),itb(k),ite(k),l,tm(l),itb(l),ite(l)
  729. 21    format(4(i9,f8.2,2i6))
  730. 20    continue
  731.     write(8,11)
  732.     do 22 i=1,nd
  733.     j=ider(3,i)
  734.     l=ipv(j)
  735. 22    write(8,24)ider(1,i),ider(2,i),(pv(m,j),m=1,l)
  736. 24    format('the derivative of',i4,' wrt',i4,' is ',15a1)
  737.     write(8,11)
  738.     write(8,40)(i,ivn(i),ivp(i),i=1,itetop)
  739. 40    format(10(1x,i5,2i3))
  740.     return
  741.     end
  742.  
  743.  
  744.     subroutine expnd(jb,jp,jans)
  745. $include: 'pfcom.for'
  746.     ia=ittb(jb)
  747.     if(ia.eq.0)call bomb(3,'expnd1  ')
  748.     if(tm(ia).ne.1.d0)call bomb(3,'expnd2  ')
  749.     ic=itb(ia)
  750.     if(ivn(ic).ne.2 .or. ivp(ic).ne.0)call bomb(3,'expnd3  ')
  751.     ic=ia+1
  752.     id=itte(jb)
  753. 5    call cpynt(ia,ia,1.d0)
  754.     if(ic.le.id)goto 10
  755. 6    call defsum(jans,nt,nt)
  756.     return
  757. 10    ig=nt
  758.     ie=ittb(jp)
  759.     tp=tm(ie)
  760.     if(tp.eq.0.d0)goto 6
  761.     jprvb=ig
  762.     jprve=ig
  763.     tprv=1.d0
  764.     do 20 kp=1,nb
  765.     if(dabs(tp).lt.1.d-11)goto 25
  766.     tprv=tprv*tp/kp
  767.     tmt(kp)=tprv
  768.     tp=tp-1.d0
  769.     itbt(kp)=nt+1
  770.     call mult(jprvb,jprve,ic,id)
  771.     itet(kp)=nt
  772.     if(itbt(kp).gt.nt)goto 25
  773.     jprvb=itbt(kp)
  774. 20    jprve=nt
  775.     kp=nb+1
  776. 25    kp=kp-1
  777.     do 40 i=1,kp
  778.     ih=itbt(i)
  779.     ii=itet(i)
  780.     tp=tmt(i)
  781.     do 30 il=ih,ii
  782. 30    tm(il)=tm(il)*tp
  783. 40    continue
  784.     call defsum(jans,ig,ii)
  785.     call order(jans)
  786.     return
  787.     end
  788.  
  789.  
  790.     subroutine extrct
  791. $include: 'pfcom.for'
  792.     call nxt(2,nchpl,' ',j1)
  793.     call srch(j1,nchpl,'-',j6)
  794.     if(j6.ne.-1)goto 300
  795.     num=1
  796.     do 1 i=2,10
  797.     if(a(j1-1).eq.hn(i))num=i-1
  798. 1    continue
  799.     call numv(j1,nchpl,j3,j4,j5)
  800.     call nxt(j4,nchpl,'=',j6)
  801.     call nxt(j6,nchpl,'m',j7)
  802.     call numv(j7+1,nchpl,j8,j9,j10)
  803.     k=ittb(j10)
  804.     if(k.eq.0)call bomb(3,'extrct1 ')
  805.     kk=tm(k)+.1
  806.     if(kk.gt.0)goto 15
  807. 12    write(8,14)tm(k),(a(i),i=j8,j9)
  808. 14    format('bad term number  ',d12.3,5x,15a1)
  809.     call bomb(3,'extrct2 ')
  810. 15    call nxt(j9+1,nchpl,'o',j11)
  811.     call numv(j11+2,nchpl,j12,j13,j14)
  812.     kl=itte(j14)+1-ittb(j14)
  813.     if(kk.gt.kl)goto 12
  814.     km=ittb(j14)+kk-1
  815.     nt1=nt+1
  816.     ktop=min0(km+num-1,itte(j14))
  817.     call cpynt(km,ktop,1.d0)
  818.     call defsum(j5,nt1,nt)
  819.     call outqu(j5)
  820.     return
  821. 300    call numv(j1,nchpl,j2,j3,j4)
  822.     call nxt(j3,j6,'=',j5)
  823.     call nxt(j5,j6,'t',j9)
  824.     call nxt(j9,j6,' ',j10)
  825.     call nxtnb(j10,j6,j11)
  826.     call srch(j11,j6,' ',j12)
  827.     if(j12.eq.-1)j12=j6
  828.     call getint(j11,j12-1,j13)
  829.     call srchnb(j6+1,nchpl,j8)
  830.     call srch(j8,nchpl,' ',j14)
  831.     call getint(j8,j14-1,j15)
  832.     call nxt(j14,nchpl,'o',j19)
  833.     if(a(j19+1).ne.'f' .and. a(j19+2).ne.' ')
  834.      +    call bomb(3,'extrct5 ')
  835. 330    call numv(j19+2,nchpl,j16,j17,j18)
  836.     if(ittb(j18).eq.0)call bomb(3,'extrct6')
  837.     if(j15.gt.itte(j18)-ittb(j18)+1)j15=itte(j18)-ittb(j18)+1
  838.     if(j15.lt.j13)call bomb(3,'extrct7 ')
  839.     nt1=nt+1
  840.     call cpynt(ittb(j18)+j13-1,ittb(j18)+j15-1,1.d0)
  841.     call defsum(j4,nt1,nt)
  842.     call outqu(j4)
  843.     return
  844.     end
  845.  
  846.  
  847.     subroutine fmtc(it,if,il,posch)
  848. $include: 'pfcom.for'
  849.     character*1 aa,pa,posch
  850.     common/pr1/pa(232),aa(22)
  851.     character*127 cdum
  852.     common/fmt/cdum
  853.     il=22
  854.     tt=tm(it)
  855.     if(tt.ne.0.)goto 3
  856.     aa(9)='0'
  857.     if=9
  858.     il=9
  859.     return
  860. 3    ab=tm(it)*64.d0
  861.     if(ab.gt.0.d0)j=ab+.5d0
  862.     if(ab.lt.0.d0)j=ab-.5d0
  863.     t=j/64.d0
  864.     ac=dmax1(1.d0,dabs(tm(it)))
  865.     if(dabs((t-tt)/ac).lt.1.d-11)tm(it)=t
  866.     tt=tm(it)
  867.     if(dabs(tt).gt.9999.9d0)goto 16
  868. 6    format(22a1)
  869. 14    if(dabs(tt).gt.1.d0)goto 20
  870.     if(dabs(tt).gt..001d0)goto 24
  871. 16    write(cdum,18)tt
  872. 18    format(d22.13)
  873.     read(cdum,6)aa
  874.     ig=22
  875. 181    ig=ig-1
  876.     if(aa(ig).ne.'E' .and. aa(ig).ne.'e')goto 181
  877.     ih=ig
  878. 182    ih=ih-1
  879.     if(aa(ih).eq.'0')goto 182
  880.     do 183 i=ig,22
  881. 183    aa(ih+1+i-ig)=aa(i)
  882.     il=ih+1+22-ig
  883.     goto 32
  884. 20    write(cdum,22)tt
  885. 22    format(f22.12)
  886.     goto 28
  887. 24    write(cdum,26)tt
  888. 26    format(f22.16)
  889. 28    read(cdum,6)aa
  890. 30    if(aa(il).ne.'0')goto 32
  891.     il=il-1
  892.     goto 30
  893. 32    do 34 if=1,il
  894.     if(aa(if).ne.' ')goto 36
  895. 34    continue
  896. 36    if(tt.gt.0.d0)if=if-1
  897.     if(tt.gt.0.d0)aa(if)=posch
  898.     return
  899.     end
  900.  
  901.  
  902.     subroutine fort(it,ib)
  903. $include: 'pfcom.for'
  904.     character*1 aa,pa
  905.     common/pr1/pa(232),aa(22)
  906.     character*127 cdum
  907.     common/fmt/cdum
  908.     if(np.eq.0)return
  909.     j=ib
  910.     call fmtc(it,if,il,'+')
  911.     do 42 i=if,il
  912.     ib=ib+1
  913. 42    a(ib)=aa(i)
  914.     ic=itb(it)
  915.     id=ite(it)
  916.     do 60 iv=ic,id
  917.     ip=iabs(ivp(iv))
  918.     in=ivn(iv)
  919.     ie=ipv(in)
  920.     iq=0
  921.     if(ip.eq.1)goto 47
  922.     write(cdum,43)ip
  923. 43    format(i20)
  924.     read(cdum,6)aa
  925. 6    format(22a1)
  926.     il=20
  927.     do 44 if=1,20
  928.     if(aa(if).ne.' ')goto 45
  929. 44    continue
  930. 45    iq=il+1-if+2
  931. 47    if(iq+ie+1.le.72)goto 52
  932.     j=0
  933.     write(8,50)(a(ich),ich=1,ib)
  934. 50    format(80a1)
  935.     do 51 i=1,80
  936. 51    a(i)=' '
  937.     a(6)='.'
  938.     ib=11
  939. 52    ib=ib+1
  940.     a(ib)='*'
  941.     if(ivp(iv).lt.0)a(ib)='/'
  942.     do 53 i=1,ie
  943.     ib=ib+1
  944. 53    a(ib)=pv(i,in)
  945.     if(iq.eq.0)goto 60
  946.     ib=ib+1
  947.     a(ib)='*'
  948.     ib=ib+1
  949.     a(ib)='*'
  950.     do 54 i=if,il
  951.     ib=ib+1
  952. 54    a(ib)=aa(i)
  953. 60    continue
  954.     k=0
  955.     if(a(j+2).eq.'1' .and. a(j+3).eq.'.' .and. a(j+4).eq.'*')k=1
  956.     if(k.eq.0)goto 61
  957.     a(j+2)=' '
  958.     a(j+3)=a(j+1)
  959.     a(j+4)=' '
  960.     a(j+1)=' '
  961. 61    write(8,50)(a(ich),ich=1,ib)
  962.     return
  963.     end
  964.  
  965.  
  966.     subroutine fortrn
  967. $include: 'pfcom.for'
  968.     call nxt(1,nchpl,' ',j1)
  969.     call numv(j1+1,nchpl,j2,j3,j4)
  970.     if(np.eq.0 .or. ittb(j4).eq.0)goto 10
  971.     do 1 i=1,nchpl
  972. 1    a(i)=' '
  973.     ia=ipv(j4)
  974.     ib=6
  975.     do 2 i=1,ia
  976.     ib=ib+1
  977. 2    a(ib)=pv(i,j4)
  978.     ib=ib+1
  979.     a(ib)='='
  980.     ic=ittb(j4)
  981.     id=itte(j4)
  982.     j=0
  983.     do 4 it=ic,id
  984.     call fort(it,ib)
  985.     do 3 i=1,nchpl
  986. 3    a(i)=' '
  987.     a(6)='.'
  988.     ib=7
  989.     j=j+1
  990.     if(j.lt.10)goto 4
  991.     a(6)=' '
  992.     call copy(pv(1,j4),a(7),ia)
  993.     a(7+ia)='='
  994.     call copy(pv(1,j4),a(8+ia),ia)
  995.     ib=7+ia+ia
  996.     j=0
  997. 4    continue
  998. 10    return
  999.     end
  1000.  
  1001.  
  1002.     subroutine getint(j1,j2,j3)
  1003. c  on return j3=integer value of a(j1)-a(j2)
  1004.     j3=0
  1005.     if(j2.lt.j1)call bomb(3,'getint1 ')
  1006.     do 1 i=j1,j2
  1007.     call intval(i,j)
  1008. 1    j3=10*j3+j
  1009.     return
  1010.     end
  1011.  
  1012.  
  1013.     subroutine getexp(j1,ie,k,if)
  1014. c  begin looking at a(j1), exponent put in ie, a(k) is next nonblank
  1015. c     (or k=-1 if only blanks are found)
  1016. $include: 'pfcom.for'
  1017.     ie=1
  1018.     call srchnb(j1,nchpl,k)
  1019.     if(k.eq.-1)return
  1020.     if(a(k).eq.'*'.and. a(k+1).ne.'*')return
  1021.     ise=1
  1022.     if(a(k).eq.'*'.and. a(k+1).eq.'*')goto 35
  1023.     if(a(k).eq.'/')return
  1024.     call srchch(j1,nchpl,j2)
  1025.     if(j2.eq.k .and. if.ne.0)return
  1026.     if(j2.ne.k)goto 9
  1027. 32    call numv(j2,nchpl,j3,j4,j5)
  1028.     if(ittb(j5).eq.0)call bomb(3,'getexp10')
  1029.     tmittb=tm(ittb(j5))
  1030.     ie=tmittb+.5d0
  1031.     if(tmittb.lt.0.d0)ie=tmittb-.5d0
  1032.     re=ie
  1033.     if(dabs(re-tmittb).gt.1.d-11)call bomb(3,'getexp11')
  1034.     k=-1
  1035.     if(j4.ge.nchpl-1)return
  1036.     call srchnb(j4+1,nchpl,k)
  1037.     return
  1038. 9    if(a(k).eq.'+')goto 20
  1039.     if(a(k).eq.'-')goto 10
  1040.     if(a(k).eq.'^')goto 30
  1041.     j3=k
  1042. 1    me=0
  1043. 5    call intval(j3,j4)
  1044.     me=10*me+j4
  1045.     j3=j3+1
  1046.     if(j3.gt.nchpl)goto 70
  1047.     if(a(j3).eq.' ')goto 50
  1048.     if(a(j3).eq.'*')goto 60
  1049.     if(a(j3).eq.'/')goto 60
  1050.     goto 5
  1051. 10    ise=-1
  1052. 20    if(k+1.gt.nchpl)call bomb(3,'getexp2 ')
  1053.     call nxtnb(k+1,nchpl,j3)
  1054.     if(j3.eq.-1)call bomb(3,'getexp4 ')
  1055.     goto 1
  1056. 30    call nxtnb(k+1,nchpl,j3)
  1057. 31    k=j3
  1058.     j2=j3
  1059.     call srchch(j3,nchpl,j4)
  1060.     if(j3.eq.j4)goto 32
  1061.     if(a(k).eq.'-')goto 10
  1062.     if(a(k).eq.'+')goto 20
  1063.     goto 1
  1064. 35    if(k.ge.nchpl-1)call bomb(3,'getexp6 ')
  1065.     call nxtnb(k+2,nchpl,j3)
  1066.     if(j3.eq.-1)call bomb(3,'getexp8 ')
  1067.     goto 31
  1068. 50    ie=me*ise
  1069.     call srchnb(j3,nchpl,k)
  1070.     return
  1071. 60    k=j3
  1072.     if(j3.eq.nchpl)call bomb(3,'getexp9 ')
  1073.     ie=me*ise
  1074.     return
  1075. 70    ie=me*ise
  1076.     k=-1
  1077.     return
  1078.     end
  1079.  
  1080.  
  1081.     function icomp(ia,ib)
  1082. $include: 'pfcom.for'
  1083. c  if ia should precede ib icomp=1
  1084. c  if ib should precede ia icomp=-1
  1085. c  if terms are the same except possibly coefficients, icomp=0
  1086.     ja=itb(ia)
  1087.     jb=itb(ib)
  1088.     ka=ite(ia)
  1089.     kb=ite(ib)
  1090. 1    if(ivn(ja)-ivn(jb))2,8,9
  1091. 2    if(ivp(ja))3,10,4
  1092. 4    icomp=-1
  1093.     return
  1094. 5    if(jb.gt.kb)goto 7
  1095. 3    icomp=1
  1096.     return
  1097. 6    ja=ja+1
  1098. 11    jb=jb+1
  1099.     if(ja.gt.ka)goto 5
  1100.     if(jb.gt.kb)goto 4
  1101.     goto 1
  1102. 7    icomp=0
  1103.     return
  1104. 8    if(ivp(ja)-ivp(jb))3,6,4
  1105. 9    if(ivp(jb))4,11,3
  1106. 10    ja=ja+1
  1107.     if(ja.gt.ka)goto 3
  1108.     goto 1
  1109.     end
  1110.  
  1111.  
  1112.     subroutine iftest(go)
  1113. $include: 'pfcom.for'
  1114.     character*1 aj1,aj2
  1115.     go=0.
  1116.     call nxt(2,nchpl,' ',j1)
  1117.     call nxt(j1,nchpl,'.',j)
  1118.     call numv(j1,j-1,j3,j4,j5)
  1119.     aj1=a(j+1)
  1120.     aj2=a(j+2)
  1121.     if(aj1.eq.'n')goto 20
  1122.     if(aj1.eq.'l' .or. aj1.eq.'g')goto 22
  1123.     if(aj1.ne.'e')call bomb(3,'iftest1 ')
  1124.     if(aj2.ne.'q')call bomb(3,'iftest2 ')
  1125. 4    if(a(j+3).ne.'.')call bomb(3,'iftest3 ')
  1126.     call numv(j+4,nchpl,j6,j7,j8)
  1127.     if(j5.ne.j8)goto 10
  1128.     if(aj1.eq.'n')return
  1129. 5    call nxtch(j7+1,nchpl,j9)
  1130.     do 6 j=j9,nchpl
  1131.     a(j+1-j9)=a(j)
  1132. 6    a(j)=' '
  1133.     go=1.d0
  1134.     call copy(a,ap,nchpl)
  1135.     return
  1136. 10    k1=itb(ittb(j5))
  1137.     k2=ite(itte(j5))
  1138.     k3=itb(ittb(j8))
  1139.     k4=ite(itte(j8))
  1140.     t5=tm(ittb(j5))
  1141.     t8=tm(ittb(j8))
  1142.     if(aj1.eq.'l')goto 23
  1143.     if(aj1.eq.'g')goto 25
  1144.     if(k2-k1.ne.k4-k3)goto 30
  1145.     k=k3
  1146.     do 11 i=k1,k2
  1147.     if(ivn(i).ne.ivn(k))goto 30
  1148.     if(ivp(i).ne.ivp(k))goto 30
  1149. 11    k=k+1
  1150.     k1=ittb(j5)
  1151.     k2=itte(j5)
  1152.     k3=ittb(j8)
  1153.     k4=itte(j8)
  1154.     if(k2-k1.ne.k4-k3)goto 30
  1155.     k=k3
  1156.     do 12 i=k1,k2
  1157.     den=dabs(tm(i))+dabs(tm(k))
  1158.     if(den.lt.1.d-11)goto 12
  1159.     if(dabs(tm(i)-tm(k))/den.gt.1.d-11)goto 30
  1160. 12    k=k+1
  1161.     if(aj1.eq.'n')return
  1162.     goto 5
  1163. 20    if(a(j+2).ne.'e')call bomb(3,'iftest4 ')
  1164.     goto 4
  1165. 22    if(aj2.eq.'e' .or. aj2.eq.'t')goto 4
  1166.     call bomb(5,'iftest5 ')
  1167. 23    if(aj2.eq.'e')goto 24
  1168.     if(t5.lt.t8)goto 5
  1169.     return
  1170. 24    if(t5.le.t8)goto 5
  1171.     return
  1172. 25    if(aj2.eq.'e')goto 26
  1173.     if(t5.gt.t8)goto 5
  1174.     return
  1175. 26    if(t5.ge.t8)goto 5
  1176.     return
  1177. 30    if(aj1.eq.'n')goto 5
  1178.     return
  1179.     end
  1180.  
  1181.  
  1182.     subroutine incnt
  1183. $include: 'pfcom.for'
  1184.     nt=nt+1
  1185.     maxnt=max0(nt,maxnt)
  1186.     if(maxnt.gt.nterm)call bomb(4,'nterm   ')
  1187.     return
  1188.     end
  1189.  
  1190.  
  1191.     subroutine intval(j1,i)
  1192. c  on return i=integer value of a(j1)
  1193. $include: 'pfcom.for'
  1194.     do 1 j=1,10
  1195.     if(a(j1).eq.hn(j))goto 2
  1196. 1    continue
  1197.     i=-1
  1198.     return
  1199. 2    i=j-1
  1200.     return
  1201.     end
  1202.  
  1203.  
  1204.     subroutine monitr
  1205. $include: 'pfcom.for'
  1206.     call readin
  1207.     ntp=0
  1208.     tprev=tnow
  1209. 10    if(ha(1).ne.'a')call bomb(3,'monitr1 ')
  1210.     if(ittb(2).ne.0)call bomb(3,'useorder')
  1211.     if(ha(1).eq.'a')goto 11
  1212.     if(ns.lt.2)goto 14
  1213.     do 20 ms=2,ns
  1214.     is=iosi(ms-1)
  1215.     js=iosi(ms)
  1216.     if(is.eq.0 .or. js.eq.0)call bomb(3,'mon1    ')
  1217.     if(ios(is).ne.ms-1 .or. ios(js).ne.ms)call bomb(3,'mon2    ')
  1218.     if(ittb(is).eq.0 .or. ittb(js).eq.0)call bomb(3,'mon3    ')
  1219.     ia=itte(is)
  1220.     ib=ittb(js)
  1221.     if(ite(ia)+1.ne.itb(ib))call bomb(3,'mon4    ')
  1222.     do 16 it=ib,itte(js)
  1223.     if(itb(it).gt.ite(it))call bomb(3,'mon5    ')
  1224.     if(dabs(tm(it)).gt.1.d-11)goto 15
  1225.     if(itb(it).ne.ite(it))call bomb(3,'mon6    ')
  1226.     if(ivn(itb(it)).ne.2)call bomb(3,'mon7    ')
  1227.     if(itb(it).gt.ite(it))call bomb(3,'mon8    ')
  1228. 15    do 17 i=itb(it),ite(it)
  1229.     if(ivn(i).eq.1)call bomb(3,'mon9    ')
  1230.     if(ittb(ivn(i)).ne.0)call bomb(3,'mon10   ')
  1231.     if(i.eq.itb(it))goto 17
  1232.     if(ivp(i).eq.0)call bomb(3,'mon11   ')
  1233.     if(ivn(i-1).ge.ivn(i))call bomb(3,'mon12   ')
  1234. 17    continue
  1235. 16    continue
  1236. 20    continue
  1237. 14    if(ns.eq.0)call bomb(3,'mon13   ')
  1238.     ic=iosi(ns)
  1239.     if(itte(ic).ne.nt)call bomb(3,'mon14   ')
  1240.     if(ite(itte(ic)).ne.itetop)call bomb(3,'mon15   ')
  1241. 11    call read(ib)
  1242.     if(ib.eq.1 .and. np.ge.1)call prntap
  1243.     if(ib.eq.1)goto 10
  1244.     call datime
  1245.     t=tnow-tprev
  1246.     ttot=tnow-tbeg
  1247.     tprev=tnow
  1248.     ntmntp=nt-ntp
  1249.     if(np.ge.3)write(8,3)t,ttot,ntmntp,nt,itetop,maxnt,maxite
  1250. 3    format(10x,'t=',f6.2,' total t=',f8.1,' ntmntp=',i5,' nt=',i5,
  1251.      +       ' itetop=',i6,' maxnt=',i5,' maxite=',i6)
  1252.     if(np.ge.2)call prntap
  1253. 5    ntp=nt
  1254.     call copy(a,ap,nchpl)
  1255.     if(a(1).eq.'c')goto 300
  1256.     if(a(1).eq.'d')goto 400
  1257.     if(a(1).eq.'e')goto 500
  1258.     if(a(1).eq.'f')goto 600
  1259.     if(a(1).eq.'g')goto 700
  1260.     if(a(1).eq.'i')goto 900
  1261.     if(a(1).eq.'l')goto 10
  1262.     if(a(1).eq.'m')goto 1300
  1263.     if(a(1).eq.'n')goto 1400
  1264.     if(a(1).eq.'o')goto 1500
  1265.     if(a(1).eq.'p')goto 1600
  1266.     if(a(1).eq.'r')goto 1800
  1267.     if(a(1).eq.'s')goto 1900
  1268.     if(a(1).eq.'t')goto 2000
  1269.     if(a(1).eq.'z')goto 2600
  1270.     call bomb(3,'bad cmnd')
  1271. 300    if(a(2).eq.'o')goto 310
  1272.     if(a(2).ne.'a')call bomb(3,'monitr2 ')
  1273.     call call
  1274.     goto 10
  1275. 310    if(a(3).eq.'m')goto 320
  1276.     if(a(3).ne.'u')call bomb(3,'monitr3 ')
  1277.     call count
  1278.     goto 10
  1279. 320    call comput
  1280.     call pack
  1281.     goto 10
  1282. 400    if(a(2).eq.'u')goto 470
  1283.     if(a(2).eq.'i')goto 320
  1284.     if(a(2).ne.'e')call bomb(3,'monitr4 ')
  1285.     if(a(3).eq.'l')goto 450
  1286.     if(a(3).eq.'f')goto 440
  1287.     call bomb(3,'monitr5 ')
  1288. 440    call dfine
  1289.     goto 10
  1290. 450    call delet
  1291.     call pack
  1292.     goto 10
  1293. 470    call dump
  1294.     goto 10
  1295. 500    if(a(2).eq.'n')return
  1296.     if(a(2).ne.'x')call bomb(3,'monitr6 ')
  1297.     call extrct
  1298.     goto 10
  1299. 600    if(a(2).eq.'l')goto 610
  1300.     call fortrn
  1301.     goto 10
  1302. 610    call readin
  1303.     goto 10
  1304. 700    call nxt(1,nchpl,' ',j9)
  1305.     call numl(j9,nchpl,j10,j11,j12)
  1306.     if(linlab(j12).eq.-1)write(8,710)
  1307.     if(linlab(j12).eq.-1)call bomb(3,'monitr7 ')
  1308. 710    format('label not found')
  1309.     line=linlab(j12)-1
  1310.     goto 10
  1311. 900    call iftest(go)
  1312.     if(go.eq.1.d0)goto 5
  1313.     goto 10
  1314. 1300    call srch(1,nchpl,'=',j)
  1315.     call srchnb(j+1,nchpl,k)
  1316.     call getexp(k,maxe,j,0)
  1317.     goto 10
  1318. 1400    call nxt(3,nchpl,'=',j1)
  1319.     call getexp(j1+1,j3,j2,0)
  1320.     if(a(2).eq.'d')goto 1405
  1321.     if(a(2).eq.'p')goto 1410
  1322.     if(a(2).eq.'b')goto 1415
  1323.     call bomb(3,'monitr8 ')
  1324. 1405    if(a(3).eq.'u')goto 1406
  1325.     call bomb(3,'monitr9 ')
  1326. 1406    ndump=j3
  1327.     goto 10
  1328. 1410    np=j3
  1329.     goto 10
  1330. 1415    nb=j3
  1331.     goto 10
  1332. 1500    if(nv.ne.1)call bomb(3,'putfirst')
  1333.     call ordset
  1334.     goto 10
  1335. 1600    if(a(3).eq.'i')goto 1620
  1336.     if(a(3).ne.'o')call bomb(3,'monitr10')
  1337. 1605    call read(ib)
  1338.     if(np.ge.2)call prntap
  1339.     if(a(1).ne.'r' .or. a(2).ne.'e' .or. a(3).ne.'t')goto 1605
  1340.     goto 10
  1341. 1620    call nxt(4,nchpl,' ',j1)
  1342.     call numv(j1,nchpl,j2,j3,j4)
  1343.     if(np.eq.1)write(8,1)
  1344.     if(np.eq.1)call prntap
  1345. 1    format(130a1)
  1346.     call print(j4)
  1347.     goto 10
  1348. 1800    if(a(2).ne.'e' .or. a(3).ne.'t')call bomb(3,'return? ')
  1349.     line=linret
  1350.     goto 10
  1351. 1900    call sub
  1352.     call pack
  1353.     goto 10
  1354. 2000    if(a(2).eq.'a')goto 2010
  1355.     call thderv
  1356.     goto 10
  1357. 2010    call taylor
  1358.     goto 10
  1359. 2600    call numv(2,nchpl,j1,j2,j3)
  1360.     call zsub(j3)
  1361.     call pack
  1362.     goto 10
  1363.     end
  1364.  
  1365.  
  1366.     subroutine mult(ia,ib,ic,id)
  1367. c  multiply terms ia-ib by terms ic-id
  1368. c  result goes to nt+1,.....
  1369. $include: 'pfcom.for'
  1370.     common/nwt/iitb,iite,jitb,jite,ntb,if
  1371.     if(ia.eq.0 .or. ib.eq.0 .or. ic.eq.0 .or. id.eq.0)
  1372.      +   call bomb(3,'mult1   ')
  1373.     ntb=nt+1
  1374.     if=0
  1375.     do 20 i=ia,ib
  1376.     iitb=itb(i)
  1377.     iite=ite(i)
  1378.     do 10 j=ic,id
  1379.     jitb=itb(j)
  1380.     jite=ite(j)
  1381.     if(ivn(iitb)+ivn(jitb).eq.4 .and. ivp(iitb)+ivp(jitb).gt.maxe)
  1382.      +                                                       goto 20
  1383.     call sub2(tm(i)*tm(j))
  1384. 10    continue
  1385. 20    continue
  1386.     return
  1387.     end
  1388.  
  1389.  
  1390.     subroutine niltrm
  1391. $include: 'pfcom.for'
  1392.     k=itetop+1
  1393.     call incnt
  1394.     tm(nt)=0.d0
  1395.     itb(nt)=k
  1396.     call nite(k)
  1397.     ivn(k)=2
  1398.     ivp(k)=0
  1399.     return
  1400.     end
  1401.  
  1402.  
  1403.     subroutine nite(j)
  1404. $include: 'pfcom.for'
  1405.     itetop=j
  1406.     ite(nt)=j
  1407.     maxite=max0(maxite,j)
  1408.     if(maxite.gt.ntot)call bomb(4,'ntot    ')
  1409.     return
  1410.     end
  1411.  
  1412.  
  1413.     subroutine numl(j1,j2,j3,j4,j5)
  1414. c  search j1-j2  find name in j3-j4, j5 is index of label
  1415. $include: 'pfcom.for'
  1416.     character*1 c
  1417.     if(j2.lt.j1)call bomb(3,'numl    ')
  1418.     call nxtch(j1,j2,j3)
  1419.     do 10 j=j3,j2
  1420.     c=a(j)
  1421.     do 5 k=1,26
  1422.     if(c.eq.ha(k) .or. c.eq.hc(k))goto 10
  1423. 5    continue
  1424.     do 6 k=1,10
  1425.     if(c.eq.hn(k))goto 10
  1426. 6    continue
  1427.     if(c.eq.' ' .or. c.eq.'+' .or. c.eq.'-')goto 12
  1428.     if(c.eq.'*' .or. c.eq.'/' .or. c.eq.'=' .or. c.eq.'^')goto 12
  1429. 10    continue
  1430.     j4=j2
  1431.     goto 14
  1432. 12    j4=j-1
  1433. 14    l=j4-j3+1
  1434.     do 18 i=1,labels
  1435.     if(ipl(i).ne.l)goto 18
  1436.     do 16 j=1,l
  1437.     if(pl(j,i).ne.a(j+j3-1))goto 18
  1438. 16    continue
  1439.     j5=i
  1440.     return
  1441. 18    continue
  1442.     labels=labels+1
  1443.     if(labels.gt.nlab)call bomb(4,'nlab    ')
  1444.     j5=labels
  1445.     ipl(j5)=l
  1446.     call copy(a(j3),pl(1,j5),l)
  1447.     return
  1448.     end
  1449.  
  1450.  
  1451.     subroutine numv(j1,j2,j3,j4,j5)
  1452. c  search j1-j2  find name in j3-j4, j5 is index of variable
  1453. $include: 'pfcom.for'
  1454.     character*1 c
  1455.     if(j2.lt.j1)call bomb(3,'numv    ')
  1456.     call nxtch(j1,j2,j3)
  1457.     do 10 j=j3,j2
  1458.     c=a(j)
  1459.     do 5 k=1,26
  1460.     if(c.eq.ha(k) .or. c.eq.hc(k))goto 10
  1461. 5    continue
  1462.     do 6 k=1,10
  1463.     if(c.eq.hn(k))goto 10
  1464. 6    continue
  1465.     if(c.eq.' ' .or. c.eq.'+' .or. c.eq.'-')goto 12
  1466.     if(c.eq.'*' .or. c.eq.'/' .or. c.eq.'=' .or. c.eq.'^')goto 12
  1467. 10    continue
  1468.     j4=j2
  1469.     goto 14
  1470. 12    j4=j-1
  1471. 14    l=j4-j3+1
  1472.     do 28 ii=ns,1,-1
  1473.     i=iosi(ii)
  1474.     if(ipv(i).ne.l)goto 28
  1475.     do 26 j=1,l
  1476.     if(pv(j,i).ne.a(j+j3-1))goto 28
  1477. 26    continue
  1478.     j5=i
  1479.     return
  1480. 28    continue
  1481.     do 18 i=1,nv
  1482.     if(ittb(i).ne.0)goto 18
  1483.     if(ipv(i).ne.l)goto 18
  1484.     do 16 j=1,l
  1485.     if(pv(j,i).ne.a(j+j3-1))goto 18
  1486. 16    continue
  1487.     j5=i
  1488.     return
  1489. 18    continue
  1490.     nv=nv+1
  1491.     j5=nv
  1492.     if(nv.gt.nvar)call bomb(4,'nvar    ')
  1493.     ipv(nv)=l
  1494.     call copy(a(j3),pv(1,nv),l)
  1495.     return
  1496.     end
  1497.  
  1498.  
  1499.     subroutine nxt(j1,j2,h,j3)
  1500.     character*1 h
  1501.     call srch(j1,j2,h,j3)
  1502.     if(j3.ne.-1)return
  1503.     write(8,1)h,j1,j2
  1504. 1    format(' nxt did not find ',a1,' in ',2i4)
  1505.     call bomb(3,'nxt     ')
  1506.     return
  1507.     end
  1508.  
  1509.  
  1510.     subroutine nxtch(j1,j2,j3)
  1511.     call srchch(j1,j2,j3)
  1512.     if(j3.eq.-1)call bomb(3,'nxtch   ')
  1513.     return
  1514.     end
  1515.  
  1516.  
  1517.     subroutine nxtnb(j1,j2,j3)
  1518.     call srchnb(j1,j2,j3)
  1519.     if(j3.eq.-1)call bomb(3,'nxtnb   ')
  1520.     end
  1521.  
  1522.  
  1523.     subroutine order(is)
  1524. $include: 'pfcom.for'
  1525.     ina=ittb(is)
  1526.     inb=itte(is)
  1527.     if(inb.lt.ina .or. ina.eq.0)call bomb(3,'order1  ')
  1528.     it=ina
  1529. 2    ia=itb(it)
  1530.     if(dabs(tm(it)).le.1.d-11)goto 41
  1531.     ib=ite(it)
  1532. 3    if(ib.eq.ia)goto 25
  1533.     i=ia
  1534. 4    if(ivp(i).eq.0)goto 20
  1535.     i=i+1
  1536.     if(i.le.ib)goto 4
  1537.     if(ia.eq.ib)goto 25
  1538. 23    i=ia+1
  1539. 5    if(ivn(i)-ivn(i-1))10,7,6
  1540. 6    i=i+1
  1541.     if(i.le.ib)goto 5
  1542.     goto 40
  1543. 7    ivp(i-1)=ivp(i-1)+ivp(i)
  1544.     if(i.eq.ib)goto 9
  1545.     ib=ib-1
  1546.     do 8 j=i,ib
  1547.     ivn(j)=ivn(j+1)
  1548. 8    ivp(j)=ivp(j+1)
  1549.     if(ivp(i-1).ne.0)goto 5
  1550.     goto 3
  1551. 9    ib=ib-1
  1552.     if(ivp(ib).ne.0)goto 40
  1553.     goto 3
  1554. 10    if(i.eq.ia+1)goto 19
  1555.     ic=i-2
  1556. 11    if(ivn(i)-ivn(ic))12,13,16
  1557. 12    if(ic.eq.ia)goto 18
  1558.     ic=ic-1
  1559.     goto 11
  1560. 13    ivp(ic)=ivp(ic)+ivp(i)
  1561.     if(i.eq.ib)goto 15
  1562.     ib=ib-1
  1563.     do 14 j=i,ib
  1564.     ivn(j)=ivn(j+1)
  1565. 14    ivp(j)=ivp(j+1)
  1566.     if(ivp(ic).eq.0)goto 3
  1567.     goto 5
  1568. 15    ib=ib-1
  1569.     if(ivp(ic).ne.0)goto 40
  1570.     goto 3
  1571. 16    ip=ivp(i)
  1572.     in=ivn(i)
  1573.     do 17 j=i,ic+2,-1
  1574.     ivn(j)=ivn(j-1)
  1575. 17    ivp(j)=ivp(j-1)
  1576.     ivn(ic+1)=in
  1577.     ivp(ic+1)=ip
  1578.     goto 6
  1579. 18    ic=ia-1
  1580.     goto 16
  1581. 19    ip=ivp(i)
  1582.     in=ivn(i)
  1583.     ivp(i)=ivp(ia)
  1584.     ivn(i)=ivn(ia)
  1585.     ivp(ia)=ip
  1586.     ivn(ia)=in
  1587.     goto 6
  1588. 20    if(i.eq.ib)goto 22
  1589.     ib=ib-1
  1590.     do 21 j=i,ib
  1591.     ivn(j)=ivn(j+1)
  1592. 21    ivp(j)=ivp(j+1)
  1593.     goto 4
  1594. 22    ib=ib-1
  1595.     if(ib.gt.ia)goto 23
  1596.     ib=ia
  1597. 25    if(ivp(ia).eq.0)ivn(ia)=2
  1598. 40    ite(it)=ib
  1599.     if(ivn(ia).ne.2 .or. ivp(ia).le.maxe)goto 52
  1600. 41    call bkup(it,inb)
  1601.     goto 53
  1602. 52    it=it+1
  1603. 53    if(it.le.inb)goto 2
  1604. c  sort terms ina-inb
  1605.     if(ina.lt.inb)goto 70
  1606.     if(ina.eq.inb)goto 82
  1607. 51    inb=ina
  1608.     tm(ina)=0.
  1609.     k=itb(ina)
  1610.     ite(ina)=k
  1611.     ivn(k)=2
  1612.     ivp(k)=0
  1613.     goto 82
  1614. 70    it=ina
  1615. 68    it=it+1
  1616. 69    if(it.gt.inb)goto 80
  1617.     if(dabs(tm(it)).le.1.d-11)goto 75
  1618.     if(icomp(ina,it))71,73,74
  1619. 71    call tzap(it,ina)
  1620.     goto 68
  1621. 73    tm(ina)=tm(ina)+tm(it)
  1622. 75    call bkup(it,inb)
  1623.     goto 69
  1624. 74    if(icomp(it-1,it))62,61,68
  1625. 61    tm(it-1)=tm(it-1)+tm(it)
  1626.     call bkup(it,inb)
  1627.     goto 69
  1628. 62    la=ina
  1629.     lb=it-1
  1630. 63    if(lb-la.eq.1)goto 64
  1631.     lc=(la+lb)/2
  1632.     if(icomp(lc,it))65,66,67
  1633. 64    call tzap(it,lb)
  1634.     goto 68
  1635. 65    lb=lc
  1636.     goto 63
  1637. 66    tm(lc)=tm(lc)+tm(it)
  1638.     call bkup(it,inb)
  1639.     goto 69
  1640. 67    la=lc
  1641.     goto 63
  1642. 80    if(ina.gt.inb)goto 51
  1643. 82    call defsum(is,ina,inb)
  1644.     return
  1645.     end
  1646.  
  1647.  
  1648.     subroutine ordset
  1649. $include: 'pfcom.for'
  1650.     call nxt(2,nchpl,' ',j1)
  1651.     call srchnb(j1,nchpl,j2)
  1652.     if(j2.eq.-1)goto 20
  1653.     if(a(j2).eq.';')goto 20
  1654.     ntsav=nt
  1655.     itesav=itetop
  1656.     call term(j1)
  1657.     nt=ntsav
  1658.     itetop=itesav
  1659.     return
  1660. 20    if=1
  1661.     do 10 il=1,300
  1662.     call read(ib)
  1663.     if(np.ge.2)call prntap
  1664.     if(ib.eq.1)return
  1665.     is=1
  1666. 2    call srchch(is,nchpl,j1)
  1667.     if(j1.eq.-1)goto 10
  1668.     call numv(j1,nchpl,j2,j3,j4)
  1669.     if=if+1
  1670.     if(j4.lt.if)call bomb(3,'ordset1 ')
  1671.     if(j4.eq.if)goto 8
  1672.     do 5 i=nv,if,-1
  1673.     ipv(i)=ipv(i-1)
  1674.     do 4 ic=1,nch
  1675. 4    pv(ic,i)=pv(ic,i-1)
  1676. 5    continue
  1677.     call copy(a(j2),pv(1,if),j3-j2+1)
  1678.     ipv(if)=j3-j2+1
  1679. 8    is=j3+1
  1680.     if(is.le.nchpl)goto 2
  1681. 10    continue
  1682.     return
  1683.     end
  1684.  
  1685.  
  1686.     subroutine outqu(j5)
  1687. $include: 'pfcom.for'
  1688. 1    do 6 it=ittb(j5),itte(j5)
  1689.     do 5 iv=itb(it),ite(it)
  1690.     if(ittb(ivn(iv)).ne.0)goto 9
  1691. 5    continue
  1692. 6    continue
  1693.     return
  1694. 9    jv=ivn(iv)
  1695.     call subs(jv,jv,j5)
  1696.     goto 1
  1697.     end
  1698.  
  1699.  
  1700.     subroutine pack
  1701. $include: 'pfcom.for'
  1702.     if(minns.le.ns)goto 9
  1703.     is=iosi(ns)
  1704.     nt=itte(is)
  1705.     itetop=ite(nt)
  1706.     minns=ns+1
  1707.     return
  1708. 9    if(minns.le.0)call bomb(3,'pack4   ')
  1709.     js=minns
  1710.     kt=0
  1711.     kv=0
  1712.     if(js.eq.1)goto 13
  1713.     kt=itte(iosi(js-1))
  1714.     kv=ite(kt)
  1715. 13    do 30 ms=js,ns
  1716.     is=iosi(ms)
  1717.     if(is.le.0)call bomb(3,'pack6   ')
  1718.     ia=ittb(is)
  1719.     ib=itte(is)
  1720.     lv=0
  1721.     lt=0
  1722.     do 14 it=ia,ib
  1723.     if(dabs(tm(it)).lt.1.d-11)goto 14
  1724.     lt=lt+1
  1725.     if(lt.gt.ntermt)call bomb(4,'ntermt  ')
  1726.     tmt(lt)=tm(it)
  1727.     itbt(lt)=lv+1
  1728.     do 12 i=itb(it),ite(it)
  1729.     if(ivp(i).ne.0)goto 31
  1730.     do 32 j=itb(it),ite(it)
  1731.     if(ivp(j).ne.0)goto 12
  1732. 32    continue
  1733.     lv=lv+1
  1734.     ifnt(lv)=2
  1735.     ifpt(lv)=0
  1736.     goto 33
  1737. 31    lv=lv+1
  1738.     ifnt(lv)=ivn(i)
  1739.     ifpt(lv)=ivp(i)
  1740. 12    continue
  1741. 33    itet(lt)=lv
  1742.     if(lv.gt.ntott)call bomb(4,'ntott   ')
  1743. 14    continue
  1744.     if(lv.ne.0)goto 15
  1745.     kt=kt+1
  1746.     ittb(is)=kt
  1747.     itte(is)=kt
  1748.     kv=kv+1
  1749.     tm(kt)=0.d0
  1750.     itb(kt)=kv
  1751.     ite(kt)=kv
  1752.     ivn(kv)=2
  1753.     ivp(kv)=0
  1754.     goto 30
  1755. 15    ittb(is)=kt+1
  1756.     do 18 it=1,lt
  1757.     kt=kt+1
  1758.     tm(kt)=tmt(it)
  1759.     itb(kt)=kv+1
  1760.     do 16 i=itbt(it),itet(it)
  1761.     kv=kv+1
  1762.     ivn(kv)=ifnt(i)
  1763. 16    ivp(kv)=ifpt(i)
  1764. 18    ite(kt)=kv
  1765.     itte(is)=kt
  1766. 30    continue
  1767.     nt=kt
  1768.     itetop=kv
  1769.     minns=ns+1
  1770.     return
  1771.     end
  1772.  
  1773.  
  1774.     subroutine print(ii)
  1775. $include: 'pfcom.for'
  1776.     common/pr1/pa(232),aa(22)
  1777.     character*1 aa,pa,bl,sm,cf
  1778.     data bl,sm,cf/' ','-','^'/
  1779.     if(np.eq.0)return
  1780.     ia=ittb(ii)
  1781.     if(ia.eq.0)write(8,2)
  1782.     if(ia.eq.0)return
  1783. 2    format('there is nothing to print')
  1784.     ib=itte(ii)
  1785.     nst=21
  1786.     do 40 it=ia,ib
  1787.     n=nst
  1788.     do 4 i=1,232
  1789. 4    pa(i)=' '
  1790.     ka=itb(it)
  1791.     kb=ite(it)
  1792.     do 30 k=ka,kb
  1793.     n=n+1
  1794.     pa(n)=bl
  1795.     if=ivn(k)
  1796.     jp=ivp(k)
  1797.     ifl=ipv(if)
  1798.     do 20 j=1,ifl
  1799.     n=n+1
  1800. 20    pa(n)=pv(j,if)
  1801.     if(jp.eq.1)goto 30
  1802.     n=n+1
  1803.     if(jp.ge.0)pa(n)=cf
  1804.     if(jp.lt.0)pa(n)=sm
  1805.     iap=iabs(jp)
  1806.     n=n+1
  1807.     j=mod(iap,10)
  1808.     pa(n)=hn(j+1)
  1809.     if(iap.le.9)goto 30
  1810.     j=iap/10
  1811.     if(j.ge.10)call bomb(3,'print2  ')
  1812.     n=n+1
  1813.     pa(n)=pa(n-1)
  1814.     pa(n-1)=hn(j+1)
  1815. 30    continue
  1816.     if(n.gt.232)call bomb(3,'print3  ')
  1817.     call fmtc(it,if,il,' ')
  1818.     if(il-if.eq.2 .and. aa(if+1).eq.'1' .and. aa(il).eq.'.')il=if
  1819.     do 32 i=if,il
  1820. 32    pa(nst-il+i)=aa(i)
  1821.     j=it+1-ia
  1822. 33    if(j.ge.1000)j=j-1000
  1823.     if(j.ge.1000)goto 33
  1824.     write(8,34)j,(pa(k),k=1,n)
  1825. 34    format(i3,2x,127a1,/,27x,105a1)
  1826. 40    continue
  1827.     return
  1828.     end
  1829.  
  1830.  
  1831.     subroutine prntap
  1832. $include: 'pfcom.for'
  1833.     do 1 n=nchpl,1,-1
  1834.     if(ap(n).ne.' ')goto 2
  1835. 1    continue
  1836.     n=1
  1837. 2    write(8,3)(ap(i),i=1,n)
  1838. 3    format(131a1)
  1839.     return
  1840.     end
  1841.  
  1842.  
  1843.     subroutine read(ib)
  1844. $include: 'pfcom.for'
  1845. c  ib=0 on return if line has nonblanks
  1846. c  ib=1 on return if line is all blanks
  1847.     line=line+1
  1848.     ib=inpba(line)
  1849.     ie=inpea(line)
  1850.     il=ie-ib+1
  1851.     do 1 i=1,nchpl
  1852. 1    a(i)=' '
  1853.     call copy(inpl(ib),a,il)
  1854.     call copy(a,ap,nchpl)
  1855.     i=1
  1856.     if(a(1).eq.'*')goto 8
  1857.     do 6 i=1,nchpl
  1858.     if(a(i).eq.';')goto 8
  1859. 6    continue
  1860.     goto 14
  1861. 8    do 10 j=i,nchpl
  1862. 10    a(j)=' '
  1863. 14    ib=0
  1864.     do 16 i=1,nchpl
  1865.     if(a(i).ne.' ')goto 18
  1866. 16    continue
  1867. 17    ib=1
  1868.     return
  1869. 18    if(i.eq.1)return
  1870.     do 20 j=i,nchpl
  1871.     a(j+1-i)=a(j)
  1872. 20    a(j)=' '
  1873.     return
  1874.     end
  1875.  
  1876.  
  1877.     subroutine readin
  1878. $include: 'pfcom.for'
  1879.     do 1 i=1,nlab
  1880. 1    linlab(i)=-1
  1881.     line=0
  1882.     lines=0
  1883.     labels=0
  1884.     iinpl=1
  1885. 5    read(7,6)a
  1886. 6    format(127a1)
  1887.     do 7 k=nchpl,1,-1
  1888.     if(a(k).ne.' ')goto 10
  1889. 7    continue
  1890.     k=1
  1891. 10    lines=lines+1
  1892.     lintot=lintot+1
  1893.     if(lines.gt.ninl)call bomb(4,'ninl    ')
  1894.     if(k.eq.1)goto 15
  1895.     call srchnb(1,k,m)
  1896.     if(a(m).eq.'p' .and. a(m+1).eq.'r' .and. a(m+2).eq.'o')goto 11
  1897.     if(a(m).ne.'l' .or. a(m+1).ne.'a' .or. a(m+2).ne.'b')goto 15
  1898. 11    call nxt(m+2,nchpl,' ',j1)
  1899.     call nxtch(j1,nchpl,j3)
  1900.     call nxt(j3,nchpl,' ',j4)
  1901.     call numl(j1,nchpl,j3,j4,j5)
  1902.     if(linlab(j5).ne.-1)write(8,12)(a(i),i=j3,j4)
  1903. 12    format('duplicate label: ',15a1)
  1904.     if(linlab(j5).ne.-1)call bomb(3,'readin1 ')
  1905.     linlab(labels)=lines
  1906. 15    inpba(lines)=iinpl
  1907.     inpea(lines)=iinpl+k-1
  1908.     if(inpea(lines).gt.ninch)call bomb(4,'ninch   ')
  1909.     call copy(a,inpl(iinpl),k)
  1910.     iinpl=iinpl+k
  1911.     if(k.eq.1)goto 5
  1912.     if(a(m).eq.'f' .and. a(m+1).eq.'l' .and. a(m+2).eq.'u')goto 16
  1913.     if(a(m).ne.'e' .or. a(m+1).ne.'n' .or. a(m+2).ne.'d')goto 5
  1914.     if(a(m+3).ne.' ' .or. a(m+4).ne.'i' .or. a(m+5).ne.'n')goto 5
  1915. 16    write(8,20)lines,labels,iinpl
  1916. 20    format('number of lines in input file=',i5,' labels='
  1917.      +      ,i2,'  characters=',i6)
  1918.     return
  1919.     end
  1920.  
  1921.  
  1922.     subroutine recip(j13,j4)
  1923. $include: 'pfcom.for'
  1924.     ia=ittb(j13)
  1925.     iz=itte(j13)
  1926.     k=itb(ia)
  1927.     if(ivn(k).eq.2 .and. ivp(k).ne.0 .and. ia.ne.iz)then
  1928.         write(8,2)
  1929.         call bomb(3,'recip1  ')
  1930. 2    format('this simple program cannot process the denominator.',
  1931.      +     '  denominators should have 1 term without variable 1.')
  1932.     endif
  1933.     ja=nt+1
  1934.     call cpynt(ia,ia,1.d0)
  1935.     tm(ja)=1.d0/tm(ja)
  1936.     do 4 i=itb(ja),ite(ja)
  1937. 4    ivp(i)=-ivp(i)
  1938.     if(iz.eq.ia)call defsum(j4,ja,nt)
  1939.     if(iz.eq.ia)return
  1940.     if(maxe.gt.2)write(8,6)
  1941.     if(maxe.gt.2)call bomb(3,'recip2  ')
  1942. 6    format('this simple program cannot divide if maxe is',
  1943.      +      ' greater than 2')
  1944.     ib=ia+1
  1945.     k=itb(ib)
  1946.     if(ivn(k).ne.2 .or. ivp(k).lt.1)write(8,2)
  1947.     if(ivn(k).ne.2 .or. ivp(k).lt.1)call bomb(3,'recip3  ')
  1948.     if(ivp(k).eq.2)goto 14
  1949.     ic=ib
  1950.     if(ib.eq.iz)goto 10
  1951. 8    k=itb(ic+1)
  1952.     if(ic.eq.iz .or. ivp(k).eq.2)goto 10
  1953.     ic=ic+1
  1954.     goto 8
  1955. 10    jb=nt+1
  1956.     call cpynt(ib,ic,1.d0)
  1957.     jc=nt
  1958.     if(ic.eq.iz)goto 12
  1959. 11    jd=nt+1
  1960.     call cpynt(ic+1,iz,1.d0)
  1961.     je=nt
  1962.     goto 16
  1963. 12    jd=nt+1
  1964.     je=jd
  1965.     itb(jd)=itetop+1
  1966.     call incnt
  1967.     call nite(itb(jd))
  1968.     k=itetop
  1969.     ivn(k)=2
  1970.     ivp(k)=2
  1971.     tm(jd)=0.d0
  1972.     goto 16
  1973. 14    jb=nt+1
  1974.     jc=jb
  1975.     itb(jb)=itetop+1
  1976.     call incnt
  1977.     call nite(itb(jb))
  1978.     k=itetop
  1979.     ivn(k)=2
  1980.     ivp(k)=1
  1981.     tm(jb)=0.d0
  1982.     ic=ib-1
  1983.     goto 11
  1984. 16    jf=nt+1
  1985.     call mult(ja,ja,jb,jc)
  1986.     jg=nt
  1987.     if(maxe.eq.1)goto 18
  1988.     jh=nt+1
  1989.     call mult(jf,jg,jf,jg)
  1990.     ji=nt
  1991.     jj=nt+1
  1992.     call mult(ja,ja,jd,je)
  1993.     jk=nt
  1994.     jl=nt+1
  1995.     call cpynt(jf,jg,-1.d0)
  1996.     call cpynt(jh,ji,1.d0)
  1997.     call cpynt(jj,jk,-1.d0)
  1998.     jm=nt
  1999.     call cpynt(ja,ja,1.d0)
  2000.     call mult(ja,ja,jl,jm)
  2001.     nt1=jm+1
  2002.     goto 24
  2003. 18    do 19 i=jf,jg
  2004. 19    tm(i)=-tm(i)
  2005.     call cpynt(ja,ja,1.d0)
  2006.     call mult(ja,ja,jf,jg)
  2007.     nt1=jg+1
  2008. 24    call defsum(j4,nt1,nt)
  2009.     call order(j4)
  2010.     return
  2011.     end
  2012.  
  2013.  
  2014.     subroutine swap(j4)
  2015. $include: 'pfcom.for'
  2016.     if=ittb(j4)
  2017.     ig=ios(j4)
  2018.     if(if.ne.0)minns=min0(minns,ios(j4))
  2019.     ios(j4)=ios(1)
  2020.     iosi(ios(j4))=j4
  2021.     ittb(j4)=ittb(1)
  2022.     itte(j4)=itte(1)
  2023.     if(if.eq.0)goto 5
  2024.     ios(1)=ig
  2025.     iosi(ig)=1
  2026.     ittb(1)=if
  2027.     itte(1)=if
  2028.     ite(if)=itb(if)
  2029.     call pack
  2030.     return
  2031. 5    call niltrm
  2032.     ittb(1)=0
  2033.     call defsum(1,nt,nt)
  2034.     return
  2035.     end
  2036.  
  2037.  
  2038.     subroutine srch(j1,j2,h,j3)
  2039. $include: 'pfcom.for'
  2040.     character*1 h
  2041.     do 1 j3=j1,j2
  2042.     if(a(j3).eq.h)return
  2043. 1    continue
  2044.     j3=-1
  2045.     return
  2046.     end
  2047.  
  2048.  
  2049.     subroutine srchch(j1,j2,j3)
  2050. $include: 'pfcom.for'
  2051.     do 2 j3=j1,j2
  2052.     do 1 i=1,26
  2053.     if(a(j3).eq.ha(i))return
  2054.     if(a(j3).eq.hc(i))return
  2055. 1    continue
  2056. 2    continue
  2057.     j3=-1
  2058.     end
  2059.  
  2060.  
  2061.     subroutine srchnb(j1,j2,j3)
  2062. $include: 'pfcom.for'
  2063.     do 1 j3=j1,j2
  2064.     if(a(j3).ne.' ')return
  2065. 1    continue
  2066.     j3=-1
  2067.     return
  2068.     end
  2069.  
  2070.  
  2071.     subroutine sub
  2072. $include: 'pfcom.for'
  2073.     character*3 aj4,aj8
  2074.     call nxt(1,nchpl,' ',j1)
  2075.     call nxtnb(j1,nchpl,j2)
  2076.     call nxtch(j1,nchpl,j3)
  2077.     if(j2.ne.j3)goto 10
  2078.     call numv(j1,nchpl,j2,j3,j4)
  2079.     aj4='sum'
  2080.     if(ittb(j4).eq.0)aj4='var'
  2081. 1    call nxt(j3+1,nchpl,'f',j5)
  2082.     if(a(j5+1).ne.'o')call bomb(3,'readin2 ')
  2083.     if(a(j5+2).ne.'r')call bomb(3,'readin4 ')
  2084.     if(a(j5+3).ne.' ')call bomb(3,'readin6 ')
  2085.     call numv(j5+3,nchpl,j6,j7,j8)
  2086.     aj8='var'
  2087.     if(ittb(j8).gt.0)aj8='sum'
  2088.     call nxt(j7+1,nchpl,'i',j9)
  2089.     if(a(j9+1).ne.'n')call bomb(3,'readin8 ')
  2090.     if(a(j9+2).ne.' ')call bomb(3,'readin10')
  2091.     call numv(j9+2,nchpl,j10,j11,j12)
  2092.     minns=min0(minns,ios(j12))
  2093.     if(aj4.eq.'num')goto 35
  2094.     if(aj4.eq.'var')goto 30
  2095.     if(aj8.eq.'sum')goto 50
  2096.     call subs(j4,j8,j12)
  2097.     return
  2098. 10    call coef(j1,j4,c)
  2099.     aj4='num'
  2100.     j3=j4-1
  2101.     goto 1
  2102. 30    if(aj8.eq.'sum')goto 40
  2103.     call subv(j4,j8,j12)
  2104.     return
  2105. 35    if(aj8.eq.'sum')call bomb(3,'readin12')
  2106.     call subn(c,j8,j12)
  2107.     return
  2108. 40    call subr(j8,j4,j12)
  2109.     return
  2110. 50    a(71)='A'
  2111.     a(72)='?'
  2112.     a(73)='0'
  2113.     call numv(71,73,j15,j16,j17)
  2114.     call subr(j8,j17,j12)
  2115.     call subs(j4,j17,j12)
  2116.     return
  2117.     end
  2118.  
  2119.  
  2120.     subroutine subn(c,j8,j12)
  2121. $include: 'pfcom.for'
  2122.     if(ittb(j12).eq.0)call bomb(3,'subn1   ')
  2123.     if(ittb(j8).ne.0)call bomb(3,'subn2   ')
  2124.     do 8 it=ittb(j12),itte(j12)
  2125.     do 5 i=itb(it),ite(it)
  2126.     if(ivn(i).ne.j8)goto 5
  2127.     tm(it)=tm(it)*c**ivp(i)
  2128.     ivp(i)=0
  2129.     if(ite(it).eq.itb(it))goto 5
  2130.     ite(it)=ite(it)-1
  2131.     if(i.gt.ite(it))goto 5
  2132.     do 1 j=i,ite(it)
  2133.     ivn(j)=ivn(j+1)
  2134. 1    ivp(j)=ivp(j+1)
  2135. 5    continue
  2136. 8    continue
  2137.     call order(j12)
  2138.     return
  2139.     end
  2140.  
  2141.  
  2142.     subroutine subr(ia,ib,ic)
  2143. c  replace pattern=first term in sum ia by variable ib in sum ic
  2144. $include: 'pfcom.for'
  2145.     if(ittb(ia).eq.0)call bomb(3,'subr1   ')
  2146.     if(ittb(ib).ne.0)call bomb(3,'subr2   ')
  2147.     nt1=nt+1
  2148.     j1=ittb(ia)
  2149.     j2=itb(j1)
  2150.     j3=ite(j1)
  2151.     j4=ittb(ic)
  2152.     j5=itte(ic)
  2153.     do 30 it=j4,j5
  2154.     j6=itb(it)
  2155.     j7=ite(it)
  2156.     minrat=9999
  2157.     do 5 j8=j2,j3
  2158.     jfp=ivn(j8)
  2159. 1    jfm=ivn(j6)
  2160.     if(jfm.gt.jfp)goto 22
  2161.     irat=ivp(j6)/ivp(j8)
  2162.     if(jfm.eq.jfp .and. irat.ge.1)goto 4
  2163.     jt=it+1-j4
  2164.     if(j2.eq.j3 .and. np.ge.2 .and. jfm.eq.jfp .and. irat.le.-1)
  2165.      +  write(8,2)jt,(pv(i,ic),i=1,nch),(pv(i,jfp),i=1,nch),
  2166.      +            (pv(i,ia),i=1,nch)
  2167. 2    format('term',i4,' of ',15a1,' has exponent of ',15a1,
  2168.      +       ' with different sign from ',15a1)
  2169.     j6=j6+1
  2170.     if(j6.le.j7)goto 1
  2171.     goto 22
  2172. 4    minrat=min0(minrat,ivp(j6)/ivp(j8))
  2173. 5    continue
  2174.     call incnt
  2175.     tm(nt)=tm(it)/tm(j1)
  2176.     itb(nt)=itetop+1
  2177.     j6=itb(it)
  2178.     k=itb(nt)
  2179.     do 20 j=j6,j7
  2180.     do 12 j8=j2,j3
  2181.     if(ivn(j8).eq.ivn(j))goto 16
  2182. 12    continue
  2183.     ivn(k)=ivn(j)
  2184.     ivp(k)=ivp(j)
  2185.     goto 20
  2186. 16    ivn(k)=ivn(j)
  2187.     ivp(k)=ivp(j)-minrat*ivp(j8)
  2188. 20    k=k+1
  2189.     ivn(k)=ib
  2190.     ivp(k)=minrat
  2191.     call nite(k)
  2192.     goto 30
  2193. 22    call cpynt(it,it,1.d0)
  2194. 30    continue
  2195.     call defsum(ic,nt1,nt)
  2196.     call order(ic)
  2197.     return
  2198.     end
  2199.  
  2200.  
  2201.     subroutine subv(j4,j8,j12)
  2202. $include: 'pfcom.for'
  2203.     if(ittb(j12).eq.0)call bomb(3,'subv1   ')
  2204.     if(ittb(j4).ne.0 .or. ittb(j8).ne.0)call bomb(3,'subv2   ')
  2205.     do 8 it=ittb(j12),itte(j12)
  2206.     do 5 i=itb(it),ite(it)
  2207.     if(ivn(i).eq.j8)ivn(i)=j4
  2208. 5    continue
  2209. 8    continue
  2210.     call order(j12)
  2211.     return
  2212.     end
  2213.  
  2214.  
  2215.     subroutine subs(in1,in2,in3)
  2216. c  substitute sum=in1 for variable=in2 in sum=in3
  2217. c  result goes in nt+1,...
  2218. $include: 'pfcom.for'
  2219.     common/nwt/iitb,iite,jitb,jite,ntb,if
  2220. c  dimension of ittbs ittes is max jp
  2221.     dimension ittbs(30),ittes(30)
  2222.     jp=1
  2223.     ia=in2
  2224.     ib=ittb(in3)
  2225.     ic=itte(in3)
  2226.     id=ittb(in1)
  2227.     ie=itte(in1)
  2228.     if(ib.eq.0 .or. id.eq.0)call bomb(3,'subs1   ')
  2229.     ntb=nt+1
  2230.     lev1=0
  2231.     do 30 it=ib,ic
  2232.     iitb=itb(it)
  2233.     iite=ite(it)
  2234.     if(ivn(iitb).eq.2)lev1=ivp(iitb)
  2235.     if(lev1.gt.maxe)goto 31
  2236.     do 5 if=iitb,iite
  2237.     if(ivn(if).ne.ia)goto 5
  2238.     ifp=ivp(if)
  2239.     lev2=0
  2240.     if(ifp.eq.1)goto 10
  2241.     if(id.eq.ie)goto 70
  2242.     if(ifp.gt.0)goto 17
  2243.     itpr=it-ib+1
  2244.     if(np.ge.2)write(8,1)(pv(i,ia),i=1,nch),
  2245.      +       (pv(i,in3),i=1,nch),itpr
  2246. 1    format(2x,15a1,' occurs to a negative power in ',15a1,
  2247.      +      '  term ',i4)
  2248. 5    continue
  2249. c  here copy term it to nt+1
  2250.     if=0
  2251.     jitb=0
  2252.     call sub2(tm(it))
  2253.     goto 30
  2254. c  here copy term it with substitution
  2255. 10    do 68 jt=id,ie
  2256.     jitb=itb(jt)
  2257.     jite=ite(jt)
  2258.     if(ivn(jitb).eq.2)lev2=ivp(jitb)
  2259.     if(lev1+lev2.gt.maxe)goto 30
  2260.     call sub2(tm(it)*tm(jt))
  2261. 68    continue
  2262.     goto 30
  2263. 70    if(ivn(itb(id)).eq.2)lev2=ivp(itb(id))*ifp
  2264.     if(lev1+lev2.gt.maxe)goto 30
  2265.     k=0
  2266.     do 71 i=itb(id),ite(id)
  2267.     k=k+1
  2268.     ifnt(k)=ivn(i)
  2269. 71    ifpt(k)=ivp(i)*ifp
  2270.     jitb=1
  2271.     jite=k
  2272.     call sub3(tm(it)*tm(id)**ifp)
  2273.     goto 30
  2274. c  here get in1**ifp
  2275. 17    jd=id
  2276.     je=ie
  2277.     ifpsav=ifp
  2278.     ntbsav=ntb
  2279.     iitbsv=iitb
  2280.     iitesv=iite
  2281.     ifsav=if
  2282.     ntbmu=nt
  2283.     itesav=itetop
  2284. c  test whether this power is in memory
  2285.     if(jp.eq.1)goto 80
  2286.     if(ifp.lt.jp)goto 27
  2287.     if(ittbs(jp).eq.0)goto 30
  2288.     if(ifp.eq.jp)goto 27
  2289.     jd=nt+1
  2290.     k=itetop
  2291.     do 86 i=ittbs(jp),ittes(jp)
  2292.     call incnt
  2293.     tm(nt)=tmt(i)
  2294.     itb(nt)=k+1
  2295.     do 84 j=itbt(i),itet(i)
  2296.     k=k+1
  2297.     ivn(k)=ifnt(j)
  2298. 84    ivp(k)=ifpt(j)
  2299.     call nite(k)
  2300. 86    continue
  2301.     je=nt
  2302.     goto 18
  2303. 80    ks=0
  2304.     nts=0
  2305. 18    jpb=nt+1
  2306.     ntbtmu=nt
  2307.     call mult(id,ie,jd,je)
  2308.     if(nt.eq.ntbtmu)then
  2309.         jp=jp+1
  2310.         ittbs(jp)=0
  2311.         nt=ntbmu
  2312.         ntb=ntbsav
  2313.         itetop=itesav
  2314.         goto 30
  2315.         endif
  2316.     jpe=nt
  2317.     jp=jp+1
  2318.     if(jp.gt.30)call bomb(4,'max jp  ')
  2319.     jd=jpb
  2320.     je=jpe
  2321.     ittbs(jp)=nts+1
  2322.     do 21 i=jpb,jpe
  2323.     nts=nts+1
  2324.     tmt(nts)=tm(i)
  2325.     itbt(nts)=ks+1
  2326.     do 20 j=itb(i),ite(i)
  2327.     ks=ks+1
  2328.     ifnt(ks)=ivn(j)
  2329. 20    ifpt(ks)=ivp(j)
  2330.     itet(nts)=ks
  2331.     if(ks.gt.ntott)call bomb(4,'ntott   ')
  2332.     if(nts.gt.ntermt)call bomb(4,'ntermt  ')
  2333. 21    continue
  2334.     ittes(jp)=nts
  2335.     if(jp.ne.ifp)goto 18
  2336.     nt=ntbmu
  2337.     itetop=itesav
  2338.     if=ifsav
  2339.     ifp=ifpsav
  2340.     iitb=iitbsv
  2341.     iite=iitesv
  2342.     ntb=ntbsav
  2343. 27    do 54 jt=ittbs(ifp),ittes(ifp)
  2344.     jitb=itbt(jt)
  2345.     jite=itet(jt)
  2346.     if(ifnt(jitb).eq.2)lev2=ifpt(jitb)
  2347.     if(lev1+lev2.gt.maxe)goto 30
  2348.     call sub3(tm(it)*tmt(jt))
  2349. 54    continue
  2350. 30    continue
  2351. 31    call defsum(in3,ntb,nt)
  2352.     return
  2353.     end
  2354.  
  2355.  
  2356.     subroutine sub2(coef)
  2357. $include: 'pfcom.for'
  2358.     common/nwt/iitb,iite,jitb,jite,ntb,if
  2359.     if(dabs(coef).lt.1.d-11)return
  2360.     call incnt
  2361.     k=itetop
  2362.     itb(nt)=k+1
  2363.     tm(nt)=coef
  2364.     ii=iitb
  2365.     if(jitb.eq.0)goto 15
  2366.     jj=jitb
  2367. 11    if(ii.gt.iite)goto 13
  2368.     if(ii.eq.if)goto 12
  2369.     if(jj.gt.jite)goto 15
  2370.     if(ivn(ii).lt.ivn(jj))goto 16
  2371.     if(ivn(ii).eq.ivn(jj))goto 57
  2372.     if(ivp(jj).eq.0)goto 18
  2373.     k=k+1
  2374.     ivn(k)=ivn(jj)
  2375.     ivp(k)=ivp(jj)
  2376. 18    jj=jj+1
  2377.     goto 11
  2378. 12    ii=ii+1
  2379.     goto 11
  2380. 13    if(jj.gt.jite)goto 60
  2381.     if(ivp(jj).eq.0)goto 17
  2382.     k=k+1
  2383.     ivn(k)=ivn(jj)
  2384.     ivp(k)=ivp(jj)
  2385. 17    jj=jj+1
  2386.     goto 13
  2387. 14    if(ii.gt.iite)goto 60
  2388. 15    if(ivp(ii).eq.0)goto 9
  2389.     k=k+1
  2390.     ivn(k)=ivn(ii)
  2391.     ivp(k)=ivp(ii)
  2392. 9    ii=ii+1
  2393.     if(ii.eq.if)goto 9
  2394.     goto 14
  2395. 16    if(ivp(ii).eq.0)goto 12
  2396.     k=k+1
  2397.     ivn(k)=ivn(ii)
  2398.     ivp(k)=ivp(ii)
  2399.     ii=ii+1
  2400.     goto 11
  2401. 57    if(ivp(ii)+ivp(jj).eq.0)goto 58
  2402.     k=k+1
  2403.     ivn(k)=ivn(ii)
  2404.     ivp(k)=ivp(ii)+ivp(jj)
  2405. 58    ii=ii+1
  2406.     jj=jj+1
  2407.     goto 11
  2408. 60    if(k.gt.itetop)goto 56
  2409.     k=itetop+1
  2410.     ivn(k)=2
  2411.     ivp(k)=0
  2412. 56    itesav=itetop
  2413.     call nite(k)
  2414.     if(nt.eq.ntb)goto 68
  2415.     if(icomp(ntb,nt))31,33,34
  2416. 31    call tzap(nt,ntb)
  2417.     goto 68
  2418. 33    tm(ntb)=tm(ntb)+tm(nt)
  2419.     nt=nt-1
  2420.     itetop=itesav
  2421.     goto 68
  2422. 34    if(icomp(nt-1,nt))62,61,68
  2423. 61    tm(nt-1)=tm(nt-1)+tm(nt)
  2424.     nt=nt-1
  2425.     itetop=itesav
  2426.     goto 68
  2427. 62    la=ntb
  2428.     lb=nt-1
  2429. 63    if(lb-la.eq.1)goto 64
  2430.     lc=(la+lb)/2
  2431.     if(icomp(lc,nt))65,66,67
  2432. 64    call tzap(nt,lb)
  2433.     goto 68
  2434. 65    lb=lc
  2435.     goto 63
  2436. 66    tm(lc)=tm(lc)+tm(nt)
  2437.     nt=nt-1
  2438.     itetop=itesav
  2439.     goto 68
  2440. 67    la=lc
  2441.     goto 63
  2442. 68    return
  2443.     end
  2444.  
  2445.  
  2446.     subroutine sub3(coef)
  2447. $include: 'pfcom.for'
  2448.     common/nwt/iitb,iite,jitb,jite,ntb,if
  2449.     if(dabs(coef).lt.1.d-11)return
  2450.     call incnt
  2451.     k=itetop
  2452.     itb(nt)=k+1
  2453.     tm(nt)=coef
  2454.     ii=iitb
  2455.     if(jitb.eq.0)goto 15
  2456.     jj=jitb
  2457. 11    if(ii.gt.iite)goto 13
  2458.     if(ii.eq.if)goto 12
  2459.     if(jj.gt.jite)goto 15
  2460.     if(ivn(ii).lt.ifnt(jj))goto 16
  2461.     if(ivn(ii).eq.ifnt(jj))goto 57
  2462.     if(ifpt(jj).eq.0)goto 18
  2463.     k=k+1
  2464.     ivn(k)=ifnt(jj)
  2465.     ivp(k)=ifpt(jj)
  2466. 18    jj=jj+1
  2467.     goto 11
  2468. 12    ii=ii+1
  2469.     goto 11
  2470. 13    if(jj.gt.jite)goto 60
  2471.     if(ifpt(jj).eq.0)goto 17
  2472.     k=k+1
  2473.     ivn(k)=ifnt(jj)
  2474.     ivp(k)=ifpt(jj)
  2475. 17    jj=jj+1
  2476.     goto 13
  2477. 14    if(ii.gt.iite)goto 60
  2478. 15    if(ivp(ii).eq.0)goto 9
  2479.     k=k+1
  2480.     ivn(k)=ivn(ii)
  2481.     ivp(k)=ivp(ii)
  2482. 9    ii=ii+1
  2483.     if(ii.eq.if)goto 9
  2484.     goto 14
  2485. 16    if(ivp(ii).eq.0)goto 12
  2486.     k=k+1
  2487.     ivn(k)=ivn(ii)
  2488.     ivp(k)=ivp(ii)
  2489.     ii=ii+1
  2490.     goto 11
  2491. 57    if(ivp(ii)+ifpt(jj).eq.0)goto 58
  2492.     k=k+1
  2493.     ivn(k)=ivn(ii)
  2494.     ivp(k)=ivp(ii)+ifpt(jj)
  2495. 58    ii=ii+1
  2496.     jj=jj+1
  2497.     goto 11
  2498. 60    if(k.gt.itetop)goto 56
  2499.     k=itetop+1
  2500.     ivn(k)=2
  2501.     ivp(k)=0
  2502. 56    itesav=itetop
  2503.     call nite(k)
  2504.     if(nt.eq.ntb)goto 68
  2505.     if(icomp(ntb,nt))31,33,34
  2506. 31    call tzap(nt,ntb)
  2507.     goto 68
  2508. 33    tm(ntb)=tm(ntb)+tm(nt)
  2509.     nt=nt-1
  2510.     itetop=itesav
  2511.     goto 68
  2512. 34    if(icomp(nt-1,nt))62,61,68
  2513. 61    tm(nt-1)=tm(nt-1)+tm(nt)
  2514.     nt=nt-1
  2515.     itetop=itesav
  2516.     goto 68
  2517. 62    la=ntb
  2518.     lb=nt-1
  2519. 63    if(lb-la.eq.1)goto 64
  2520.     lc=(la+lb)/2
  2521.     if(icomp(lc,nt))65,66,67
  2522. 64    call tzap(nt,lb)
  2523.     goto 68
  2524. 65    lb=lc
  2525.     goto 63
  2526. 66    tm(lc)=tm(lc)+tm(nt)
  2527.     nt=nt-1
  2528.     itetop=itesav
  2529.     goto 68
  2530. 67    la=lc
  2531.     goto 63
  2532. 68    return
  2533.     end
  2534.  
  2535.  
  2536.     subroutine switch(j,k,is)
  2537. $include: 'pfcom.for'
  2538.     do 2 it=ittb(is),itte(is)
  2539.     do 1 i=itb(it),ite(it)
  2540.     if(ivn(i).eq.j)ivn(i)=k
  2541. 1    continue
  2542. 2    continue
  2543.     return
  2544.     end
  2545.  
  2546.  
  2547.     subroutine taylor
  2548. $include: 'pfcom.for'
  2549.     character*4 star
  2550.     dimension ltay(mitay)
  2551.     data iorder,ifirst/1,0/
  2552.     if(idflg.ne.0)call dertab
  2553.     call nxt(2,nchpl,' ',j1)
  2554.     call srch(j1,nchpl,'=',j2)
  2555.     j4=j1
  2556.     if(j2.eq.-1)goto 5
  2557.     call getexp(j2+1,iorder,j3,0)
  2558.     if(iorder.lt.1)call bomb(3,' taylor1')
  2559.     return
  2560. 5    call nxt(2,nchpl,' ',j13)
  2561.     star='no'
  2562.     if(a(j13-1).eq.'*')star='yes'
  2563. 10    call nxt(j4,nchpl,'i',j3)
  2564.     if(j3.eq.-1)call bomb(3,' taylor2')
  2565.     if(a(j3-1).eq.' ' .and. a(j3+1).eq.'n' .and. a(j3+2).eq.' ')
  2566.      +       goto 12
  2567.     j4=j3+1
  2568.     if(j4.lt.nchpl)goto 10
  2569.     call bomb(3,' taylor3')
  2570. 12    do 14 j5=j3-1,j1,-1
  2571.     if(a(j5).ne.' ')goto 16
  2572. 14    continue
  2573.     call bomb(3,' taylor4')
  2574. 16    do 18 j6=j5,j1,-1
  2575.     if(a(j6).eq.' ')goto 20
  2576. 18    continue
  2577.     call bomb(3,' taylor5')
  2578. 20    call numv(j6,j5,j7,j8,j9)
  2579.     if(ittb(j9).eq.0)call bomb(3,' taylor6')
  2580.     itesav=itetop
  2581.     call term(j3+2)
  2582.     itay=itetop-itesav
  2583.     if(itay.gt.mitay)call bomb(4,'mitay   ')
  2584.     do 21 i=1,itay
  2585. 21    ltay(i)=ivn(itesav+i)
  2586.     nt=nt-1
  2587.     itetop=itesav
  2588.     call niltrm
  2589.     call defsum(1,nt,nt)
  2590.     if(ifirst.ne.0)goto 22
  2591.     a(1)='A'
  2592.     a(2)='?'
  2593.     a(3)='0'
  2594.     call numv(1,3,j14,j15,j10)
  2595.     a(3)='1'
  2596.     call numv(1,3,j14,j15,j11)
  2597.     a(3)='2'
  2598.     call numv(1,3,j14,j15,j12)
  2599. 22    ifirst=1
  2600.     nv1=nv+1
  2601.     do 34 it=ittb(j9),itte(j9)
  2602.     call cpynt(it,it,1.d0)
  2603.     call defsum(j10,nt,nt)
  2604.     do 32 ivv=1,itay
  2605.     iv=ltay(ivv)
  2606.     if(star.eq.'yes')call switch(iv,nv1,j10)
  2607.     call copys(j10,j11)
  2608.     call subn(0.d0,iv,j10)
  2609.     f=1.d0
  2610.     do 30 io=1,iorder
  2611.     call deriv(j11,iv,j12)
  2612.     call outqu(j12)
  2613.     if(tm(ittb(j12)).eq.0.d0)goto 31
  2614.     call copys(j12,j11)
  2615.     call subn(0.d0,iv,j12)
  2616.     f=f*io
  2617.     call incnt
  2618.     k=itetop+1
  2619.     tm(nt)=1.d0/f
  2620.     itb(nt)=k
  2621.     ite(nt)=k
  2622.     ivn(k)=iv
  2623.     ivp(k)=io
  2624.     call nite(k)
  2625.     nt1=nt+1
  2626.     call mult(nt,nt,ittb(j12),itte(j12))
  2627.     call defsum(j12,nt1,nt)
  2628.     call pack
  2629. 30    call add(j12,j10,1.d0,j10)
  2630. 31    if(star.eq.'yes')call switch(nv1,iv,j10)
  2631. 32    continue
  2632. 34    call add(j10,1,1.d0,1)
  2633.     call swap(j9)
  2634.     call order(j9)
  2635.     call delsum(j10)
  2636.     return
  2637.     end
  2638.  
  2639.  
  2640.     subroutine term(ist)
  2641. c  each term is on a separate line, free format with coefficient first
  2642. $include: 'pfcom.for'
  2643.     call incnt
  2644.     call coef(ist,j1,c)
  2645.     tm(nt)=c
  2646.     itb(nt)=itetop+1
  2647.     if(j1.eq.-1)then
  2648.         k=itb(nt)
  2649.         ivn(k)=2
  2650.         ivp(k)=0
  2651.         goto 10
  2652.         endif
  2653.     k=itetop
  2654. 5    ise=1
  2655.     if(a(j1).eq.'*')j1=j1+1
  2656.     if(a(j1).eq.'/')ise=-1
  2657.     if(a(j1).eq.'/')j1=j1+1
  2658.     call numv(j1,nchpl,j2,j3,j4)
  2659.     ie=1
  2660.     j5=-1
  2661.     if(j3.eq.nchpl-1 .and. a(nchpl).ne.' ')call bomb(3,'term1   ')
  2662.     if(j3.le.nchpl-2)call getexp(j3+1,ie,j5,1)
  2663.     k=k+1
  2664.     ivn(k)=j4
  2665.     ivp(k)=ie*ise
  2666.     j1=j5
  2667.     if(j5.ne.-1)goto 5
  2668. 10    call nite(k)
  2669.     return
  2670.     end
  2671.  
  2672.  
  2673.     subroutine thderv
  2674. $include: 'pfcom.for'
  2675.     idflg=1
  2676.     call nxt(2,nchpl,'o',j1)
  2677.     call nxt(j1,nchpl,' ',j2)
  2678.     call numv(j2,nchpl,j3,j4,j5)
  2679.     call nxt(j4+1,nchpl,'w',j6)
  2680.     call nxt(j6,nchpl,' ',j7)
  2681.     call numv(j7,nchpl,j8,j9,j10)
  2682.     call srch(j9+1,nchpl,'=',j11)
  2683.     if(j11.ne.-1)goto 1
  2684.     call srch(j9+1,nchpl,'i',j12)
  2685.     if(j12.eq.-1)call bomb(3,'thderv1 ')
  2686.     call nxt(j12,nchpl,' ',j11)
  2687. 1    call srchnb(j11+1,nchpl,j12)
  2688.     if(j12.eq.-1)call bomb(3,'thderv2 ')
  2689.     ig=0
  2690.     if(a(j12).eq.'0')ig=1
  2691.     if(ig.eq.1)goto 8
  2692.     call numv(j11+1,nchpl,j12,j13,j14)
  2693. 8    if=0
  2694.     do 10 i=1,nd
  2695.     if(ider(1,i).eq.j5 .and. ider(2,i).eq.j10)if=1
  2696.     if(ider(1,i).eq.j5 .and. ider(2,i).eq.j10)in=i
  2697. 10    continue
  2698.     if(if.eq.0 .and. ig.eq.1)goto 50
  2699.     if(if.eq.0)goto 14
  2700.     ja=ider(3,in)
  2701.     jb=ipv(ja)
  2702.     if(np.ge.2)write(8,11)(pv(i,ja),i=1,jb)
  2703. 11    format('REDEFINITION: previously the  derivative was = ',15a1)
  2704.     if(ig.eq.0)goto 18
  2705.     nd=nd-1
  2706.     if(in.gt.nd)goto 50
  2707.     do 28 ih=in,nd
  2708.     ider(1,ih)=ider(1,ih+1)
  2709.     ider(2,ih)=ider(2,ih+1)
  2710. 28    ider(3,ih)=ider(3,ih+1)
  2711.     goto 50
  2712. 14    nd=nd+1
  2713.     in=nd
  2714.     if(nd.gt.nder)call bomb(4,'nder    ')
  2715.     ider(1,in)=j5
  2716.     ider(2,in)=j10
  2717. 18    ider(3,in)=j14
  2718. 50    return
  2719.     end
  2720.  
  2721.  
  2722.     subroutine tzap(ka,kb)
  2723. c  zap term at ka to position kb
  2724. $include: 'pfcom.for'
  2725.     if(ka.eq.kb)return
  2726.     tmka=tm(ka)
  2727.     itbka=itb(ka)
  2728.     iteka=ite(ka)
  2729.     if(ka.gt.kb)goto 5
  2730.     do 1 i=ka,kb-1
  2731.     tm(i)=tm(i+1)
  2732.     itb(i)=itb(i+1)
  2733. 1    ite(i)=ite(i+1)
  2734.     goto 8
  2735. 5    do 6 i=ka,kb+1,-1
  2736.     tm(i)=tm(i-1)
  2737.     itb(i)=itb(i-1)
  2738. 6    ite(i)=ite(i-1)
  2739. 8    tm(kb)=tmka
  2740.     itb(kb)=itbka
  2741.     ite(kb)=iteka
  2742.     return
  2743.     end
  2744.  
  2745.  
  2746.  
  2747.     subroutine zsub(j3)
  2748.     return
  2749.     end
  2750.